5.2 Implementation
5.2.1 Plain MPC Algorithm
10 15
x
y
Figure 5.5: Stable (solid line) and unstable equilibrium (dashed line) of the uncon-trolled catalytic rod model (5.5).
display fragments of the source code which are meaningful. A complete example can be found in Appendix A.
The implementation of theMPC-POD algorithm is presented in Section 5.2.2. Here, the aim is to explain the interaction between the components rather than investi-gating the source code in detail. We complete the chapter with the presentation of the adaptive horizon MPC algorithms.
5.2.1 Plain MPC Algorithm
The principle structure of the implementation is displayed in Figure 5.6. The ab-stract classes are written in black while the example subclass is red. On the bottom of the pyramid we have the abstract classModel. It contains information about the problem setting like PDE and cost functional. The class has to provide predictions for state and adjoint equation as well as gradient information and possibly control constraints for the optimal control solver. Since we have to discretize the PDE for the simulation, this class generally requires an external ODE solver for time inte-gration and a method to discretize the spatial variable, e.g. a FEM class. However, in the case of discrete time systems or if the spatial variable is directly discretized by finite differences these classes are not necessary. An example subclass is given by CatalyticRod.
On the next level we have the abstract class Optimize. This is the centerpiece of the implementation: For the current state y(n) it solves the optimal control problem and provides the next statey(n+ 1) as well as the optimal control sequence u⋆. For this purpose it requires the information from theModelclass. An example subclass is given by BFGS which is an implementation of Algorithm 4.9.
5 Numerical Implementation
Figure 5.6: Program structure.
On the top of the pyramid we have the main program mpcpde.cpp. Here the ac-tual MPC loop takes place. Moreover, the concrete instances of the abstract classes Modeland Optimize are created and memory for the variables is allocated.
Class Model
Now, we present the abstract class Model and the concrete realization ReactDiff in detail. As already said the class provides simulations for state and adjoint equation.
Furthermore, the optimizer requires information how the gradient looks like. In Table 5.2 the abstract methods of class Modelare displayed, i.e., these methods are obligatory for the derived subclass.
The methodpredictStatesimulates the state equation with initial statexfrom initial timet tot+h for a given controlu. The solution at timet+his stored inx. More-over, the simulation for the adjoint equation is done by the method predictAdjoint.
The simulation runs from t to t+h and the solution is again stored in x. Since the adjoint equation requires information about the current state, the solution of the state equation is provided by the variable rpar. Note that we introduce a time transformation to make sure that the equation runs forward in time. Therefore, the
5.2 Implementation given state has to run backwards in time what is taken into account in the internal implementation. The last obligatory method is computeGradient which computes the current gradient by the adjoint p and the control u. The result is stored in gradient.
In the next step we want to demonstrate the structure of the presented method
Function Description variables
predictState Simulate state equation t: Initial time x: Current state u: Current control h: Simulation time predictAdjoint Simulate adjoint equation t: Initial time
x: Current adjoint rpar: Current state h: Simulation time computeGradient Compute the gradient p: Current adjoint
u: Current control
gradient: Current gradient
lambda: Regularization parameter N: Horizon
Table 5.2: Methods of the abstract classModelwhich are obligatory for every derived subclass.
for a concrete example. We investigate the derived subclass ReactDiff which is a prototype class for general reaction diffusion equations. As a concrete example we look at the one dimensional Schl¨ogl model with distributed control (5.2). In order to provide a simulation for the state equation, we have to discretize the PDE. For this easy example we use a finite difference method for the spatial variable. The method stateEquation of the class ReactDiff provides the semidiscretized PDE for the ODE solver. In Listing 5.1 the implementation of the semidiscretization is dis-played. The central finite differences are visible in the for loop. For the first and the last point the homogeneous Dirichlet condition has to take into account. The method semiFunction computes the nonlinearity of this semilinear PDE. Moreover, we observe that we have a distributed control because u acts in each discretization point independently.
Listing 5.1: Semidiscretized state equation
void R e a c t D i f f :: s t a t e E q u a t i o n (int *n , double *t , double *x , double * dx , double *u , int * ipar ) {
int dim = _ d i m e n s i o n _ s t a t e ;
double hx = LENGTH /double( dim +1);
5 Numerical Implementation
dx [0]=( - 2.0* x [0] + x [1])/ pow ( hx ,2) + s e m i F u n c t i o n ( x [0]) + u [0];
for(int i =1; i < dim -1; i ++) {
dx [ i ]=( x [i -1] -2.0* x [ i ]+ x [ i +1])/ pow ( hx ,2) + s e m i F u n c t i o n ( x [ i ]) + u [ i ];
}
dx [ dim -1]=( x [ dim -2] -2.0* x [ dim -1] )/ pow ( hx ,2) + s e m i F u n c t i o n ( x [ dim -1]) + u [ dim -1];
}
double R e a c t D i f f :: s e m i F u n c t i o n (double x ) {
f = MU *( x - pow (x ,3));
return f ; }
Again, we want to mention that we simplified the example to clarify the structure.
In the next step we present the semidiscretization of the adjoint equation. In Listing 5.2 the method adjointEquation is displayed. It can be seen that we use again cen-tral finite differences. Furthermore, we have the discrete representation of f′(y)p, where f′(y) is computed by the methodderiFunction. The last term originates from Jy =y−yd. In our example we look at yd ≡0.
Listing 5.2: Semidiscretized adjoint equation
void R e a c t D i f f :: a d j o i n t E q u a t i o n (int *n , double *t , double *x ,
double * dx , double * rpar , int * ipar ) {
int dim = _ d i m e n s i o n _ s t a t e ;
double hx = LENGTH /double( dim +1);
dx [0]=( -2.0* x [0] + x [1])/ pow ( hx ,2)
+ x [0]* d e r i F u n c t i o n ( rpar [0])+ rpar [0]
for(int i =1; i < dim -1; i ++) {
dx [ i ]=( x [i -1] -2.0* x [ i ]+ x [ i +1])/ pow ( hx ,2) + x [ i ]* d e r i F u n c t i o n ( rpar [ i ])+ rpar [ i ];
}
dx [ dim -1]=( x [ dim -2] -2.0* x [ dim -1] )/ pow ( hx ,2)
+ x [ dim -1]* d e r i F u n c t i o n ( rpar [ dim -1])+ rpar [ dim -1];
}
double R e a c t D i f f :: d e r i F u n c t i o n (double x ) {
f = MU *(1. -3.* pow (x ,2));
return f ;
5.2 Implementation
}
We want to point out that the functions described so far are only auxiliary methods for the simulation of the PDEs. The last step is to solve the semidiscretized PDE with an appropriate ODE solver. This is done by the obligatory methodspredictState and predictAdjoint, see Listing 5.3. First, the used ODE solver has to be declared in the constructor of ReactDiff (not shown in the listing) before it can be initialized by init(t,x). Afterwards the calculation is done with calc(t+h,u).
Listing 5.3: Obligatory methods predictState and predictAdjoint
void R e a c t D i f f :: p r e d i c t S t a t e ( double t , double *x , double *u , double h ) {
_odesolver - > init (t , x );
_odesolver - > calc ( t +h , u );
}
void R e a c t D i f f :: p r e d i c t A d j o i n t ( double t , double *x , double * rpar , double h ) {
_ o d e s o l v e a djo int - > init (t , x );
_ o d e s o l v e a djo int - > calc ( t +h , rpar );
}
The remaining obligatory method iscomputeGradient, where the current gradient is determined. In our exampleReactDiffthe gradient is given byJ′ =p(x, t)+λu(x, t), see Section 2.2. The resulting discretized version is displayed in Listing 5.4. We compute the gradient by controluand adjointpwhile the result is stored ingradient.
Listing 5.4: Obligatory method computeGradient
void R e a c t D i f f :: c o m p u t e G r a d i e n t ( double *p , double *u ,
double * gradient , double _lambda , int _ h o r i z o n) {
for( int j =0;j < _ h o r i z o n; j ++) {
for( int i =0; i < _ c t r l _ d i m e n s i o n ; i ++) {
g r a d i e n t[ j * _ c t r l _ d i m e n s i o n + i ]
=( p [ j * _ c t r l _ d i m e n s i o n + i ]+ _ l a m b d a* u [ j * _ c t r l _ d i m e n s i o n + i ]);
} } }
After the presentation of the obligatory methods we want to pay attention to the control constrained case. Bounds on the control can be incorporated by setting the variables control lb and control ub to the lower and upper bound. If the bounds are not defined in the derived subclass, the variables have the default values −1019 and 1019 respectively.
Finally, we want to consider the special case where we use the Newton-CG method as optimization algorithm. In Section 4.1.2 we discussed that the Newton method
5 Numerical Implementation
requires second order information. This information can also be included in the subclassReactDiff. By defining linearizedP DE =trueone can define the methods predictLinear where the linearized state equation is computed and secondDerivative where the second derivative of the semilinear term is given. Otherwise the relevant functions are determined by numerical differentiation.
The constructor of a Modelsubclass requires the state (DIMENSION X) and the control dimension (DIMENSION U). In Listing 5.5 the generation of an instance is displayed.
Listing 5.5: Initialization of a class ReactDiff.
Model * model = new R e a c t D i f f ( DIMENSION_U , D I M E N S I O N _ X );
The Modelsubclasses for the PDE examples presented in Section 5.1 can be found in Table 5.3.
Class Description
ReactDiff One dimensional Schl¨ogl model with distributed control (5.2) SemiLinBC One dimensional Schl¨ogl model with boundary control (5.3) MODELFEM FEM Schl¨ogl model with boundary control (5.4)
CatalyticRod Catalytic rod model (5.5)
Table 5.3: Example subclasses from the superclass Model.
Class Optimize
In this section we introduce the abstract class Optimize which is the superclass for our optimization algorithms. Again, we only present the public methods which are relevant for the user. In Table 5.4 these methods are displayed. The most impor-tant function is given by calc(x,u). The input variables are the current state x and an initial guess for the optimal control problem u, e.g., the shifted optimal control sequence from the previous MPC step. After its computation, the state at the next time step is stored in x while u contains the current control sequence. If the meth-ods setTolerance and setMaxIterations are not specified then the algorithm uses the default values tol = 10−6 and maxsteps = 500. The current values of the stage cost and the objective function are especially of interest for adaptive horizon MPC.
The same holds for the method resizeHorizon, where the optimization horizon can be changed in each MPC step.
In Listing 5.6 we see the initialization of a subclass BFGS. The description of the required variables is displayed in Table 5.5. If the optimization horizon changes during the MPC algorithm, the variable N has to be set to the maximum horizon.
Listing 5.6: Initialization of the class BFGS.
O P T I M I Z E * o p t i m i z e=new BFGS ( model , yd , N , lambda , T );
5.2 Implementation
Function Description variables
calc Calculation of the optimal control x: Current state u: Current control setTolerance Set prescribed accuracy tol: Tolerance
setMaxIterations Set maximum number of iterations maxsteps: Maximum steps getStageCost Return stage cost
getObjective Return objective function value
resizeHorizon Resize optimization horizon N: New horizon Table 5.4: Public methods of the abstract class Optimize.
Type Name Description Model model Example model double* yd Desired state
int N Optimization horizon double lambda Regularization parameter double T Sampling time
Table 5.5: Required variables for the constructor of a derived subclass.
In Table 5.6 we see the currently implemented subclasses of the abstract class OP-TIMIZE. The details concerning the optimization algorithms can be found in Section 4.1 while the corresponding numerical results are presented in Chapter 6. In Section
Class Description
PGM Projected Gradient Method
NCG Nonlinear Conjugate Gradient Method
BFGS BFGS Method
BFGSINV Inverse BFGS Method BFGSMF Meshfree BFGS Method NEWTONCG Newton-CG method
Table 5.6: Example subclasses from the superclass Optimize.
6.2 we distinguish between the algorithms BFGSINV I and BFGSINV II. Both are variants of BFGSINV with different initial approximations for the inverse Hessian, see Section 6.2 for details. By using the method setHessianone can choose the unit matrix (false) or the approximation from the previous MPC step (true). The default value is true because this variant turned out to be the faster one.
5 Numerical Implementation