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 $x$ collecting the vertical positions and accelerations of the masses, while the output $y$ 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, $\Delta y$ 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);

$\Delta y$ : 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 $\Delta r$

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 $\Delta r(N)$ as a function of N, that is the higher-lower sample times ratio, and $\Delta r(N)/N$ 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 $\Delta r(N)$ and $\Delta r(N)/N$. 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 $\Delta r(N)/N$

    [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 $\Delta y$ for allowing the determination of a MOAS).