Reachability Analysis and Control Synthesis for Uncertain Linear Systems in MPT

Reference:

This features were presented in: M. Kvasnica – B. Takács – J. Holaza – D. Ingole: Reachability Analysis and Control Synthesis for Uncertain Linear modeltems in MPT. In Proceedings of the 8th IFAC Symposium on Robust Control Design, Elsevier, Bratislava, Slovak Republic, 8, pp.302–307, doi.org/10.1016/j.ifacol.2015.09.474, 2015.

Uncertain linear time-invariant (LTI) system

The uncertain linear time-invariant system can be constructed as:

model = ULTISystem('A',A, 'B',B, 'E',E)

Note, the vertices of the polytopic uncertainty can be specified as a cell array, i.e., A = {A1, A2, ..., An} and B = {B1, B2, ..., Bm}. The E matrix can be omitted, in which case E = I.

Box constraints (intervals)

If simple box constraints are required, these can be set via:

model.x.min = xmin
model.x.max = xmax
model.u.min = umin
model.u.max = umax
model.d.min = dmin
model.d.max = dmax

If the constraints are not set, the default values are:

model.x.min = -Inf
model.x.max =  Inf
model.u.min = -Inf
model.u.max =  Inf
model.d.min =  0
model.d.max =  0

Polytopic constraints

If polytopic constraints are desired, they can be set by constructing a polyhedron S and running the following sequence of commands:

S = Polyhedron(H,h)
model.d.with('setConstraint')
model.d.setConstraint = S

It is important to notice that the min/max bounds can be used in combination with the polytopic constraints. In such a case, the final constraint is the intersection of the two sets. If only polytopic boundaries are set, then the final set is calculated from the given polytope and default min/max values (see Section Box constraints). Note, if only polytopic disturbance bounds D are set, the intersection with the default min/max values would be zero. In this case, we need to specify also the min/max values, e.g. as:

model.d.min = -Inf
model.d.max = Inf

Robust reachable sets

In this toolbox, the backward and forward reachable sets can be computed by performing the Pontryagin difference, Minkowski sum, and direct/inverse linear mapping of sets. Fortunately, MPT comes with a powerful polytopic library that implements all these operations and exposes them using a simple user interface.

Low-level implementation

Given polytopes P and Q. The Minkowski sum can be obtained by calling:

M = P+Q

Similarly, the Pontryagin difference is obtained by:

D = P - Q

The inverse linear mapping is obtained by:

Q = P*M 

Finally, the direct affine transformation of the polytope P is performed by:

Q = M*P 

The one-step backward reachable set of the set X is obtained using the previous operations:

preX = ((X - (E*D)) + (-B*U))*A

Note, the sets X and D have to be specified as instances of the Polyhedron class. Similarly, the one-step reachable set of the set X is obtained via:

reachX = (A*X) + (E*D)

If polytopic uncertainty is considered, then the pre-set of X is given as the intersection of pre-sets computed for the vertices of the polytopic uncertainty as follows:

preX = (X - (E*D))*A{1}
    for i = 2 : numel(A)
    preX = preX & (X - (E*D))*A{i}
end

We recall that the polytopic uncertainty is represented by the cell array of vertices A = {A1, ..., An}.

Similarly, if polytopic uncertainty is considered, the one-step reachable set does not have to be convex and is given as the union of sets since all possible values of matrix A need to be considered. An inner approximation can be obtained by considering the vertices of polytopic uncertainty:

reachX = [];
for i = 1 : numel(A)
    for j = 1 : numel(B)
        reachX = [reachX, (A{i}*X) + (B{j}*U) + (E*D)]
    end
end

User interface

MPT provides a single-command interface for the computation of one-step robust reachable sets. The type of the reachable set, i.e., forward or backward direction, needs to be specified first:

direction = 'forward'  % one-step reachable set
direction = 'backward' % one-step backward reachable set
S = model.reachableSet('direction', direction, 'X', X, 'U', U, 'D', D)

The X, U, and D arguments allow the users to specify particular sets X, U, and/or D described in the previous subsection. If one of the X, U, and/or D inputs is omitted, the toolbox uses the default values described above, see Section Box constraints.

Demo:

Here we state the full demo example:

model = ULTISystem('A', [0.5, 0; 1, -0.5])
model.d.min = [-1; -1]
model.d.max = [1, 1]
X = Polyhedron('lb', [-10;-10], 'ub', [10;10])
preX = model.reachableSet('X', X, 'direction', 'backward') % one-step backward reachable set

D = Polyhedron('lb',[-1;-1], 'ub', [1;1])
preX = model.reachableSet('X', X, 'D', D, 'direction', 'backward') % one-step backward reachable set with constrained disturbance

reachX = model.reachableSet('X', X, 'direction', 'forward') % one-step reachable set

Robust invariant set

Two types of robust invariant sets can be computed by MPT: (maximal) robust positive invariant sets for autonomous systems and (maximal) robust control invariant sets for non-autonomous uncertain systems.

The maximal robust positive/control invariant set Ω utilizes the one-step backward reachable set and is calculated by recursion. The procedure is initialized by the set X and is terminated when Ω(k+1) = Ω(k):

Omega = X
while true
    pre = model.reachableSet('X', Omega, 'direction', 'backward')
    Omega_next = pre & Omega
    if Omega_next == Omega
        break
    else
        Omega = Omega_next
    end
end

This recursive procedure is performed by the invariantSet method which automatically calculates maximal robust positive/control invariant sets for uncertain LTI systems:

Omega = model.invariantSet()

If the system described by the model object has control inputs, the method returns the maximal control invariant set. For autonomous system, the method provides the maximal positive invariant set.

The toolbox uses the default values of sets X, U and D described above, see Section Box constraints. If required, custom choices of sets can be provided using the longer syntax:

model.invariantSet('X', X, 'U', U, 'D', D, 'maxIterations', maxiter)

Demo:

Here we state the whole demo example to compute the maximal robust control invariant set:

A = [1.5, 0; 1, -1.5]
B = [1; 0]
model = ULTISystem('A', A, 'B', B)
model.x.min = [-10; -10]
model.x.max = [10; 10]
model.u.min = -5
model.u.max = 5
model.d.min = [-0.1; -0.1]
model.d.max = [0.1; 0.1]
Cinf = model.invariantSet()

Demo:

As a second example, consider an autonomous uncertain LTI system with polytopic uncertainty represented by the vertices of matrix A stated in cell array Aunc. The maximal positive invariant set is computed:

Aunc = {[0.50;1-0.5], [10;1-0.5]}
model = ULTISystem('A',Aunc)
model.x.min = [-10; -10]
model.x.max = [10; 10]
Oinf = model.invariantSet()

Robust model predictive control

In Multi-Parametric Toolbox, it is also possible to design a robust model predictive controller. The design procedure starts by specification of the prediction model as an instance of the ULTISystem. Subsequently, state/input/disturbance bounds can be provided. Next, the penalty matrices Qx and Qu can be set when specifying the form of the cost function:

model.x.penalty = QuadFunction(Qx) 
model.u.penalty = QuadFunction(Qu)    % Quadratic cost function

model.x.penalty = OneNormFunction(Qx)
model.u.penalty = OneNormFunction(Qu) % One-norm cost function

model.x.penalty = InfNormFunction(Qx)
model.u.penalty = InfNormFunction(Qu) % Infinity-norm cost function

No terminal set and zero terminal penalty are assumed by default. If required, they can be enabled as follows:

model.x.with('terminalPenalty')
model.x.terminalPenalty = QuadFunction(QN)
model.x.with('terminalSet')
model.x.terminalSet = T

Once the prediction model along with constraints, penalties and prediction horizon N is specified, a robust MPC controller can be obtained by:

ctrl = MPCController(model,N)

To obtain the optimal control action for some particular initial condition x0, use the evaluate method:

u0 = ctrl.evaluate(x0)

If required, an extended solution can be obtained by:

[u0,feasible,mpc] = ctrl.evaluate(x0)

which also returns the feasibility flag and a structure containing information about the optimal solution, such as the value of the cost function in mpc.cost and the open-loop profile of control inputs in mpc.U.

Robust MPC design DEMO:

%% DEMO
% ULTI System
A = [1.5, 0; 1, -1.5];
B = [1; 0];
model = ULTISystem('A', A, 'B', B)
model.x.min = [ -10; -10 ];
model.x.max = [  10;  10 ];
model.u.min = -5;
model.u.max =  5;
model.d.min = [ -0.1; -0.1 ];
model.d.max = [  0.1;  0.1 ];
Cinf = model.invariantSet();
% Penalty functions
Q = [ 1 0; 0 1 ];
R = [ 1 ];
model.x.penalty = QuadFunction(Q); 
model.u.penalty = QuadFunction(R);
% Terminal penalty
P = model.LQRPenalty;
model.x.with('terminalPenalty');
model.x.terminalPenalty = P;
% Terminal set
Tset = model.LQRSet;
model.x.with('terminalSet');
model.x.terminalSet = Tset;
% Prediction horizon
N = 5; 
% Robust MPC construction
robustMPC = MPCController(model,N)
% Initial condition
x0 = [-5; -3];
% Robust MPC evaluation
[u0, feasible, mpc] = robustMPC.evaluate(x0)