Example file for class HiMPC
In this demo we use the HiMPC toolbox to show how to build a decentralized hierarchical controller for linear systems subject to linear constraints on input and output variables.
Contents
Description
For this demo we use the following system consting of four balls and seven springs.
The process is composed by four masses moving vertically, each one connected by a spring to a fixed ceiling, subject to damping due to viscous friction with the environment, and connected to its neighbor mass by another spring. An input force u [Nm] can be applied to each mass by the lower-level controller. The state of the system is the vector collecting the vertical positions and accelerations of the masses, while the output
is the vector that collects the positions. The constrained masses-springs system is an extension of the example presented [1]. However in this demo both layers have a decentralized structure, in particular four controllers at the higher level independently process the user provided references and command their respective actuators.
The final objective of this class functions is to compute the maximum reference variation, element-wise, that the higher controller can apply in order to guarantee the system stability while fulfilling the constraints. Such goal is achieved by ensuring that the system state vector starting in a Maximal Output Admissible Set (MOAS) will be in a new MOAS at the next execution of the higher level controller. The MOAS is obtained by tightening the output constraints by a fixed value, that is a tuning knob of the approach. The computation of maximum reference variation is an optimization problem, which can be conveniently casted as a Mixed Integer Linear Program (MILP), and solved via suitable solvers (we strongly suggest IBM, former ILOG, Cplex).
The simulations are performed using at the higher layer a set of independent hierarchical MPC (HiMPC) which adopt a linear formulation. The upper-layer controllers must embed constraints on the generated references, to ensure stability and constraint satisfaction. In this example we impose constraints on plant outputs, namely the positions of the balls depicted in figure.
[1] D. Barcelli, A. Bemporad and G. Ripaccioli, Hierarchical MPC controller design for LTI systems, 49th Control and Decision Conference, Atlanta, Georgia, USA, 2010.
Simulation Example
close all clear variables clear globals clc
Plant model
Sample time
TL=.25;
Creates the closed-loop model of the plant with the lower layer controller
[sys N]=buildModel(TL); [A B C D] = ssdata(sys);
Number of states and inputs respectively
[Nx Ny]=size(B);
: Tightening of the output admissible set. Tuning knob.
DeltaY=.3;
Output constraints
Ycon.min = -0.3*ones(N,1); Ycon.max = 1*ones(N,1);
State constraints The output constraints have to be expressed in state coordinates in order to computer the MOAS. To avoid numerical problems the state constraints polytope is bounded, hence the speed is limited to 10 [m/s]
Xcon.min=[]; for i=1:N,Xcon.min=[Xcon.min;-DeltaY;-10];end Xcon.max=[]; for i=1:N,Xcon.max=[Xcon.max;DeltaY;10];end
coupledCons(1).H=[]; coupledCons(1).K=[]; coupledCons(2).H=[]; coupledCons(2).K=[]; coupledCons(3).H=[]; coupledCons(3).K=[]; coupledCons(4).H=[]; coupledCons(4).K=[];
The decentralization of the model is defined as:
dec(1).x=[1 2]; dec(1).y=1; dec(1).u=1; dec(1).applied=1; dec(2).x=[3 4]; dec(2).y=2; dec(2).u=2; dec(2).applied=2; dec(3).x=[5 6]; dec(3).y=3; dec(3).u=3; dec(3).applied=3; dec(4).x=[7 8]; dec(4).y=4; dec(4).u=4; dec(4).applied=4;
Create the instance of the object
himpc_obj = HiMPC(sys,dec,Ycon,DeltaY,Xcon,coupledCons);
looking for available solvers... MPT toolbox 2.6.2 initialized... Copyright (C) 2003-2006 by M. Kvasnica, P. Grieder and M. Baotic Send bug reports, questions or comments to mpt@control.ee.ethz.ch For news, visit the MPT web page at http://control.ee.ethz.ch/~mpt/ LP solver: CDD Criss-Cross QP solver: quadprog MILP solver: CPLEX9 MIQP solver: CPLEX9 Vertex enumeration: CDD Run 'mpt_studio' to start the GUI. Run 'mpt_setup' to set global parameters.
Compute the MOAS. For this purpose each sub-model is considered independent and subject to an unknown but bounded disturbance. The disturbance models the influence of the unmodeled dynamics, i.e. the state components that are not present in that sub-system. The uncertainty polytope is computed in an exact manner by considering the range of each unmodeled state component and its influence with respect to state components which belong to the current sub-model.
himpc_obj=himpc_obj.computeMOAS();
computing the MOAS... polytope object: 1-by-1 computing the MOAS... polytope object: 1-by-1 computing the MOAS... polytope object: 1-by-1 computing the MOAS... polytope object: 1-by-1
Plot both the MOAS with (blue) and without (red) the influence of the other sub-models for each sub-model
himpc_obj.plotMOAS();
nDec = 4

Solve the MILP for values of sample-time ratio from 1 to 50 and, accordingly to the ratio, compute the maximal allowable
himpc_obj=himpc_obj.computeDeltaR(50);
Time needed for the 1-th problem of the 1 subsystem Elapsed time is 0.489564 seconds. DeltaR=0.00014355 Time needed for the 2-th problem of the 1 subsystem Elapsed time is 0.331708 seconds. DeltaR=0.00058961 Time needed for the 3-th problem of the 1 subsystem Elapsed time is 0.311138 seconds. DeltaR=0.0014143 Time needed for the 4-th problem of the 1 subsystem Elapsed time is 0.324457 seconds. DeltaR=0.0026698 Time needed for the 5-th problem of the 1 subsystem Elapsed time is 0.339021 seconds. DeltaR=0.0044055 Time needed for the 6-th problem of the 1 subsystem Elapsed time is 0.321295 seconds. DeltaR=0.0067362 Time needed for the 7-th problem of the 1 subsystem Elapsed time is 0.323451 seconds. DeltaR=0.0098341 Time needed for the 8-th problem of the 1 subsystem Elapsed time is 0.328178 seconds. DeltaR=0.013972 Time needed for the 9-th problem of the 1 subsystem Elapsed time is 0.322997 seconds. DeltaR=0.019592 Time needed for the 10-th problem of the 1 subsystem Elapsed time is 0.335033 seconds. DeltaR=0.027454 Time needed for the 11-th problem of the 1 subsystem Elapsed time is 0.324962 seconds. DeltaR=0.038958 Time needed for the 12-th problem of the 1 subsystem Elapsed time is 0.302777 seconds. DeltaR=0.056989 Time needed for the 13-th problem of the 1 subsystem Elapsed time is 0.309729 seconds. DeltaR=0.088584 Time needed for the 14-th problem of the 1 subsystem Elapsed time is 0.311622 seconds. DeltaR=0.15641 Time needed for the 15-th problem of the 1 subsystem Elapsed time is 0.317342 seconds. DeltaR=0.44626 Time needed for the 16-th problem of the 1 subsystem Elapsed time is 0.303903 seconds. DeltaR=0.46066 Time needed for the 17-th problem of the 1 subsystem Elapsed time is 0.302103 seconds. DeltaR=0.46513 Time needed for the 18-th problem of the 1 subsystem Elapsed time is 0.314564 seconds. DeltaR=0.46762 Time needed for the 19-th problem of the 1 subsystem Elapsed time is 0.303286 seconds. DeltaR=0.46789 Time needed for the 20-th problem of the 1 subsystem Elapsed time is 0.308135 seconds. DeltaR=0.48174 Time needed for the 21-th problem of the 1 subsystem Elapsed time is 0.326288 seconds. DeltaR=0.51276 Time needed for the 22-th problem of the 1 subsystem Elapsed time is 0.295425 seconds. DeltaR=0.55511 Time needed for the 23-th problem of the 1 subsystem Elapsed time is 0.281844 seconds. DeltaR=0 Time needed for the 1-th problem of the 2 subsystem Elapsed time is 0.304160 seconds. DeltaR=0.00014871 Time needed for the 2-th problem of the 2 subsystem Elapsed time is 0.305246 seconds. DeltaR=0.0009689 Time needed for the 3-th problem of the 2 subsystem Elapsed time is 0.313652 seconds. DeltaR=0.0025531 Time needed for the 4-th problem of the 2 subsystem Elapsed time is 0.332407 seconds. DeltaR=0.0049828 Time needed for the 5-th problem of the 2 subsystem Elapsed time is 0.368794 seconds. DeltaR=0.0083968 Time needed for the 6-th problem of the 2 subsystem Elapsed time is 0.309247 seconds. DeltaR=0.013037 Time needed for the 7-th problem of the 2 subsystem Elapsed time is 0.305010 seconds. DeltaR=0.019292 Time needed for the 8-th problem of the 2 subsystem Elapsed time is 0.352467 seconds. DeltaR=0.027795 Time needed for the 9-th problem of the 2 subsystem Elapsed time is 0.300156 seconds. DeltaR=0.039616 Time needed for the 10-th problem of the 2 subsystem Elapsed time is 0.296006 seconds. DeltaR=0.05667 Time needed for the 11-th problem of the 2 subsystem Elapsed time is 0.363703 seconds. DeltaR=0.08274 Time needed for the 12-th problem of the 2 subsystem Elapsed time is 0.297349 seconds. DeltaR=0.12691 Time needed for the 13-th problem of the 2 subsystem Elapsed time is 0.304651 seconds. DeltaR=0.21246 Time needed for the 14-th problem of the 2 subsystem Elapsed time is 0.314012 seconds. DeltaR=0.41186 Time needed for the 15-th problem of the 2 subsystem Elapsed time is 0.314832 seconds. DeltaR=0.42408 Time needed for the 16-th problem of the 2 subsystem Elapsed time is 0.326213 seconds. DeltaR=0.43372 Time needed for the 17-th problem of the 2 subsystem Elapsed time is 0.344971 seconds. DeltaR=0.44061 Time needed for the 18-th problem of the 2 subsystem Elapsed time is 0.366107 seconds. DeltaR=0.44436 Time needed for the 19-th problem of the 2 subsystem Elapsed time is 0.326523 seconds. DeltaR=0.44427 Time needed for the 20-th problem of the 2 subsystem Elapsed time is 0.351041 seconds. DeltaR=0.46053 Time needed for the 21-th problem of the 2 subsystem Elapsed time is 0.294103 seconds. DeltaR=0.49195 Time needed for the 22-th problem of the 2 subsystem Elapsed time is 0.341920 seconds. DeltaR=0.62992 Time needed for the 23-th problem of the 2 subsystem Elapsed time is 0.297418 seconds. DeltaR=0.58686 Time needed for the 24-th problem of the 2 subsystem Elapsed time is 0.297711 seconds. DeltaR=0.64616 Time needed for the 25-th problem of the 2 subsystem Elapsed time is 0.398261 seconds. DeltaR=0.64672 Time needed for the 26-th problem of the 2 subsystem Elapsed time is 0.330789 seconds. DeltaR=0.66025 Time needed for the 27-th problem of the 2 subsystem Elapsed time is 0.277271 seconds. DeltaR=0 Time needed for the 1-th problem of the 3 subsystem Elapsed time is 0.345701 seconds. DeltaR=0.00014507 Time needed for the 2-th problem of the 3 subsystem Elapsed time is 0.351044 seconds. DeltaR=0.00094878 Time needed for the 3-th problem of the 3 subsystem Elapsed time is 0.343906 seconds. DeltaR=0.0025231 Time needed for the 4-th problem of the 3 subsystem Elapsed time is 0.381569 seconds. DeltaR=0.0049392 Time needed for the 5-th problem of the 3 subsystem Elapsed time is 0.355585 seconds. DeltaR=0.0083396 Time needed for the 6-th problem of the 3 subsystem Elapsed time is 0.339756 seconds. DeltaR=0.012966 Time needed for the 7-th problem of the 3 subsystem Elapsed time is 0.341222 seconds. DeltaR=0.019207 Time needed for the 8-th problem of the 3 subsystem Elapsed time is 0.342832 seconds. DeltaR=0.027693 Time needed for the 9-th problem of the 3 subsystem Elapsed time is 0.334668 seconds. DeltaR=0.039489 Time needed for the 10-th problem of the 3 subsystem Elapsed time is 0.333471 seconds. DeltaR=0.056503 Time needed for the 11-th problem of the 3 subsystem Elapsed time is 0.348931 seconds. DeltaR=0.082516 Time needed for the 12-th problem of the 3 subsystem Elapsed time is 0.330955 seconds. DeltaR=0.12613 Time needed for the 13-th problem of the 3 subsystem Elapsed time is 0.354350 seconds. DeltaR=0.21196 Time needed for the 14-th problem of the 3 subsystem Elapsed time is 0.386079 seconds. DeltaR=0.41165 Time needed for the 15-th problem of the 3 subsystem Elapsed time is 0.354907 seconds. DeltaR=0.4238 Time needed for the 16-th problem of the 3 subsystem Elapsed time is 0.358468 seconds. DeltaR=0.43338 Time needed for the 17-th problem of the 3 subsystem Elapsed time is 0.406512 seconds. DeltaR=0.4402 Time needed for the 18-th problem of the 3 subsystem Elapsed time is 0.422057 seconds. DeltaR=0.44387 Time needed for the 19-th problem of the 3 subsystem Elapsed time is 0.343331 seconds. DeltaR=0.44407 Time needed for the 20-th problem of the 3 subsystem Elapsed time is 0.415703 seconds. DeltaR=0.46039 Time needed for the 21-th problem of the 3 subsystem Elapsed time is 0.338005 seconds. DeltaR=0.49186 Time needed for the 22-th problem of the 3 subsystem Elapsed time is 0.341536 seconds. DeltaR=0.62993 Time needed for the 23-th problem of the 3 subsystem Elapsed time is 0.358738 seconds. DeltaR=0.58685 Time needed for the 24-th problem of the 3 subsystem Elapsed time is 0.335945 seconds. DeltaR=0.64617 Time needed for the 25-th problem of the 3 subsystem Elapsed time is 0.293247 seconds. DeltaR=0.64671 Time needed for the 26-th problem of the 3 subsystem Elapsed time is 0.342526 seconds. DeltaR=0.66024 Time needed for the 27-th problem of the 3 subsystem Elapsed time is 0.304910 seconds. DeltaR=0 Time needed for the 1-th problem of the 4 subsystem Elapsed time is 0.365161 seconds. DeltaR=0.00014858 Time needed for the 2-th problem of the 4 subsystem Elapsed time is 0.353070 seconds. DeltaR=0.00059477 Time needed for the 3-th problem of the 4 subsystem Elapsed time is 0.346131 seconds. DeltaR=0.0014344 Time needed for the 4-th problem of the 4 subsystem Elapsed time is 0.348853 seconds. DeltaR=0.0027004 Time needed for the 5-th problem of the 4 subsystem Elapsed time is 0.340505 seconds. DeltaR=0.0044576 Time needed for the 6-th problem of the 4 subsystem Elapsed time is 0.341137 seconds. DeltaR=0.0068162 Time needed for the 7-th problem of the 4 subsystem Elapsed time is 0.352919 seconds. DeltaR=0.0099506 Time needed for the 8-th problem of the 4 subsystem Elapsed time is 0.368722 seconds. DeltaR=0.014137 Time needed for the 9-th problem of the 4 subsystem Elapsed time is 0.373565 seconds. DeltaR=0.019822 Time needed for the 10-th problem of the 4 subsystem Elapsed time is 0.352014 seconds. DeltaR=0.027774 Time needed for the 11-th problem of the 4 subsystem Elapsed time is 0.340619 seconds. DeltaR=0.039408 Time needed for the 12-th problem of the 4 subsystem Elapsed time is 0.333469 seconds. DeltaR=0.057667 Time needed for the 13-th problem of the 4 subsystem Elapsed time is 0.337505 seconds. DeltaR=0.089615 Time needed for the 14-th problem of the 4 subsystem Elapsed time is 0.342310 seconds. DeltaR=0.15812 Time needed for the 15-th problem of the 4 subsystem Elapsed time is 0.380050 seconds. DeltaR=0.44992 Time needed for the 16-th problem of the 4 subsystem Elapsed time is 0.360884 seconds. DeltaR=0.46124 Time needed for the 17-th problem of the 4 subsystem Elapsed time is 0.348215 seconds. DeltaR=0.46577 Time needed for the 18-th problem of the 4 subsystem Elapsed time is 0.344696 seconds. DeltaR=0.4683 Time needed for the 19-th problem of the 4 subsystem Elapsed time is 0.349870 seconds. DeltaR=0.46854 Time needed for the 20-th problem of the 4 subsystem Elapsed time is 0.355387 seconds. DeltaR=0.48217 Time needed for the 21-th problem of the 4 subsystem Elapsed time is 0.399113 seconds. DeltaR=0.51304 Time needed for the 22-th problem of the 4 subsystem Elapsed time is 0.334506 seconds. DeltaR=0.55529 Time needed for the 23-th problem of the 4 subsystem Elapsed time is 0.303241 seconds. DeltaR=0
Plot both as a function of N, that is the higher-lower sample times ratio, and
for each sub-model
himpc_obj.plotDeltaR();


Pre-simulation considerations
The first and the forth sub-models MOAS are completely equivalent as well as the corresponding and
. The same consideration can be extended to the second and third sub-models. This is motivated by the internal symmetry of both the plant and the decentralization.
Higher level controller setup
In this section we define the higher layer controllers. We propose a model predictive control (MPC) design strategy for such a layer. The function lincon from the hybrid toolbox is used to perform that task.
warning off nDec=length(dec); for i=1:nDec
High level sample time is chosen guaranteeing the maximal allowable slope for the references, as the maximum of the ratio
[sup ind] = max(himpc_obj.DrN{i}); drOpt=himpc_obj.Dr{i}(ind);
compute the Higher level sample-time
TH(i) = ind*TL;
Sub-models definition
Ai=A(dec(i).x,dec(i).x); Bi=B(dec(i).x,dec(i).u); Ci=C(dec(i).y,dec(i).x); Di=D(dec(i).y,dec(i).u); [Ai Bi Ci Di]=ssdata(d2d(ss(Ai,Bi,Ci,Di,TL),TH(i))); nx=length(dec(i).x); nu=length(dec(i).u); ny=length(dec(i).y);
State and reference polytopes (MOAS) to be imposed in order to guarantee stability
Ko=himpc_obj.Ko; Ho=himpc_obj.Ho; Ky=himpc_obj.Ky; Hy=himpc_obj.Hy; ncx(i)=length(Ko{i}); ncu(i)=length(Ky{i}); Ci=[Ci;Ho{i};zeros(ncu(i),nx)]; Di=[Di;-Ho{i}*Ci(1:ny,:)';Hy{i}];
Controller tuning
sysOld=ss(Ai,Bi,Ci,Di,TL);
sysNew=d2d(sysOld,TH(i));
model=sysNew;
type='track';
MPC setup, 'help lincon' for more details
cost.Q = .01*eye(nx); cost.R = .01*eye(nu); cost.S=blkdiag(1*eye(ny),.001*eye(ncx(i)+ncu(i))); cost.T=.01*eye(nu); cost.rho=1e6; interval.N=2; interval.Nu=1; limits.ymin=[Ycon.min(dec(i).y);-inf*ones(ncx(i)+ncu(i),1)];
Tighted output constraint computation
KK=Ky{i}-DeltaY; limits.ymax=[Ycon.max(dec(i).y); Ko{i} ; KK]; limits.dumin=-drOpt*ones(nu,1); limits.dumax=drOpt*ones(nu,1);
Controller object creation
L{i}=lincon(model,type,cost,interval,limits,'cplex',1);
Hybrid Toolbox v.1.1.12 [February 24, 2009] - (C) 2003-2009 by Alberto Bemporad <bemporad@dii.unisi.it>
Hybrid Toolbox v.1.1.12 [February 24, 2009] - (C) 2003-2009 by Alberto Bemporad <bemporad@dii.unisi.it>
Hybrid Toolbox v.1.1.12 [February 24, 2009] - (C) 2003-2009 by Alberto Bemporad <bemporad@dii.unisi.it>
Hybrid Toolbox v.1.1.12 [February 24, 2009] - (C) 2003-2009 by Alberto Bemporad <bemporad@dii.unisi.it>
end
Simulation
Simulation time
Tsim = 60; t=0:TL:Tsim-TL;
Number of iteration of the discrete time loop
TT=Tsim/TL;
Define references
r=[ .65*ones(1,TT/2-15) .2*ones(1,TT/2+15); .8*ones(1,TT/4+10) .7*ones(1,TT*2/4-10) .4*ones(1,TT*1/4); 0.0*ones(1,TT/2) 0.7*ones(1,TT/2); 0.0*ones(1,TT)];
Initial condition
x0=[0 0 0 0 0 0 0 0]';
Initial output
u0=zeros(4,1); x=[x0 x0]; u=u0;
Simulation loop
for k=2:TT-1
For each controller
for i=1:nDec
Determine which states are used by the current controller
xx=x(dec(i).x,k);
%yy=[xx;himpc_obj.Ho{i}*xx;himpc_obj.Hy{i}*u(dec(i).u,k-1)];
Determine the references that applies to the current sub-system
rr=[r(dec(i).y,k);zeros(ncx(i)+ncu(i),1)];
If it is not time to change the reference (input given by the current controller) then use the previous one
u(dec(i).u,k) = u(dec(i).u,k-1);
Check if it is time for the higher level controller to execute
if mod(k,TH(i)/TL)==2
The controller type is 'track' the MPC will return the input increment
u(dec(i).u,k) = u(dec(i).u,k)+eval(L{i},xx,rr,u(dec(i).u,k-1));
end
end
Update the state
x(:,k+1)=A*x(:,k)+B*u(:,k);
Assume the next input to be the same as current one. That is done for plot reasons which request x and u to have the same length
u(:,k+1)=u(:,k);
end
Plot the first (continuous blue) and second (dashed blue) mass positions against their references (continuous red and dash red), respectively
figure(4); plot(t,x(1,:)','b',t,x(3,:)','--b',t,r(1,:)','r',t,r(2,:)','--r'); title('First and second mass positions')

Plot first and second inputs against their references
figure(5); plot(t,u(1,:)','b',t,u(2,:)','--b',t,r(1,:)','r',t,r(2,:)','--r') title('First and second inputs')

Plot showing all the mass trajectories against their references
figure(6); plot(t,x(1:2:7,:)','b',t,r,'r') title('All mass positions')

Plot showing all inputs (applied references) against their references
figure(7); plot(t,u','b',t,r,'r') title('All inputs')

Considerations
Figures 4 and 6 show that the hard constraint that bounds the mass upper position holds for both masses. Moreover the tracking is effective with limited overshoot. Figure 4 also show that the system is actually coupled since both masses 3 and 4, whose reference is stationary after time instant 20, have some fluctuation due to the interaction with neighbors.
Figures 5 and 7 show the references actually applied to the lower level controller. It is interesting to see how limits imposed have the effect of limiting the overshoot on the mass position. Moreover, it should be observed that the reference on the second mass (dashed blue line) does not reach the user defined reference (red continuous line) which is not in the admissible set for references (we recall that the admissible reference set is the output set tighter by the tuning know factor for allowing the determination of a MOAS).