Main /

MPT3 allows to modify the cost function and the constraints of arbitrary MPC setups via YALMIP. The functionality is implemented in the `MPCController/toYALMIP()`

and `MPCController/fromYALMIP()`

methods.

As an illustrative example, consider the following scenario:

- prediction model {$ x^+ = Ax + Bu $} with {$ A=0.9 $} and {$ B=1 $}
- state constraints {$ -5 \le x \le 5 $}
- input constraints {$ -1 \le u \le 1 $}
- unity quadratic penalties on states and inputs
- prediction horizon {$ N = 3 $}

Such an MPC problem is formulated by:

model = LTISystem('A', 0.9, 'B', 1);

model.x.min = -5;

model.x.max = 5;

model.u.min = -1;

model.u.max = 1;

model.x.penalty = QuadFunction(1);

model.u.penalty = QuadFunction(1);

N = 3;

mpc = MPCController(model, N);

model.x.min = -5;

model.x.max = 5;

model.u.min = -1;

model.u.max = 1;

model.x.penalty = QuadFunction(1);

model.u.penalty = QuadFunction(1);

N = 3;

mpc = MPCController(model, N);

Now suppose that we want the input constraints to be time-varying:

- {$ -0.5 \le u_0 \le 0.5 $} for {$u_0$} (the first element of the predicted input sequence)
- {$ -0.8 \le u_1 \le 0.8 $} for {$u_1$} (the second element)

This can be easily achieved by modifying the MPC problem via YALMIP. As a first step, you need to obtain the YALMIP representation of the MPC controller:

Y = mpc.toYALMIP();

The variable `Y`

is a structure that contains the constraints, objective, and variables.

Now we can directly modify the constraints as desired:

Y.constraints = Y.constraints + [ -0.5 <= Y.variables.u(:, 1) <= 0.5 ];

Y.constraints = Y.constraints + [ -0.8 <= Y.variables.u(:, 2) <= 0.8 ];

Y.constraints = Y.constraints + [ -0.8 <= Y.variables.u(:, 2) <= 0.8 ];

Finally, we tell the MPC controller to use the new constraints:

mpc.fromYALMIP(Y);

As an another example we show how to include a non-convex terminal set into MPC setups. The assumption here is that the non-convex set can be decomposed into a finite number of polyhedra, say `P = [P1, P2, P3, P4]`

. To tell MPT3 to use `P`

as the terminal set, use

Y = mpc.toYALMIP();

Y.constraints = Y.constraints + [ ismember( Y.variables.x(:, end), P) ];

mpc.fromYALMIP(Y);

Y.constraints = Y.constraints + [ ismember( Y.variables.x(:, end), P) ];

mpc.fromYALMIP(Y);

Here, `Y.variables.x(:, end)`

denotes the final predicted state, while the `ismember()`

command tells MPT/YALMIP to ensure that the terminal state resides in one of the elements of the polyhedron array `P`

.

Note that nonconvex set constraints lead to mixed-integer problems which are difficult to solve.

To add the constraint {$ F x_k + G u_k \le b $} for each step of the prediction horizon, do:

Y = mpc.toYALMIP();

for k = 1:size(Y.variables.u, 2)

Y.constraints = Y.constraints + [ F*Y.variables.x(:, k) + G*Y.variables.u(:, k) <= b ];

end

mpc.fromYALMIP(Y);

for k = 1:size(Y.variables.u, 2)

Y.constraints = Y.constraints + [ F*Y.variables.x(:, k) + G*Y.variables.u(:, k) <= b ];

end

mpc.fromYALMIP(Y);

If these modifications or additional constraints are not satisfactory to customize your MPC problem, it is still possible to formulate the problem completely in YALMIP. For more details on formulating the problems in YALMIP, see MPC examples in YALMIP. In the example below, it is shown how to formulate MPC problem where the upper/lower bounds on inputs are parametrized.

% parametrization of upper/lower bounds on the input using YALMIP

% dynamics x(k+1) = A*x(k) + Bu(k)

A = 0.9;

B = 1;

% penalties

Q = 5;

R = 0.1;

horizon = 5;

% variables

x = sdpvar(1,horizon+1);

u = sdpvar(1,horizon);

ub = sdpvar(1,1);

lb = sdpvar(1,1);

con = [];

obj = 0;

for i=1:horizon

% assing domain

con = [con; x(i+1) == A*x(i) + B*u(i) ];

% input constraints

con = [con; lb <= u(i) <= ub ];

% constraints on parameters

con = [con; -2 <= lb <= 0; lb <= ub; 0 <= ub <= 2];

% state constraints

con = [con; -10 <= x(i) <= 10 ];

obj = obj + norm(Q*(x(i)-1),1) + norm(R*u(i),1);

end

% construct the optimization problem with [x0; lb; ub] as parameters and u

% as the decision variables

problem = Opt(con, obj, [x(1); lb; ub], u);

% solve the problem explicitly

res = problem.solve

% evaluate the primal solution u for a value of parameter [0;-1;1]

res.xopt.feval([0;-1;1],'primal')

% evaluate the primal solution u for a value of parameter [0;-1;0.5]

res.xopt.feval([0;-1;0.5],'primal')

% dynamics x(k+1) = A*x(k) + Bu(k)

A = 0.9;

B = 1;

% penalties

Q = 5;

R = 0.1;

horizon = 5;

% variables

x = sdpvar(1,horizon+1);

u = sdpvar(1,horizon);

ub = sdpvar(1,1);

lb = sdpvar(1,1);

con = [];

obj = 0;

for i=1:horizon

% assing domain

con = [con; x(i+1) == A*x(i) + B*u(i) ];

% input constraints

con = [con; lb <= u(i) <= ub ];

% constraints on parameters

con = [con; -2 <= lb <= 0; lb <= ub; 0 <= ub <= 2];

% state constraints

con = [con; -10 <= x(i) <= 10 ];

obj = obj + norm(Q*(x(i)-1),1) + norm(R*u(i),1);

end

% construct the optimization problem with [x0; lb; ub] as parameters and u

% as the decision variables

problem = Opt(con, obj, [x(1); lb; ub], u);

% solve the problem explicitly

res = problem.solve

% evaluate the primal solution u for a value of parameter [0;-1;1]

res.xopt.feval([0;-1;1],'primal')

% evaluate the primal solution u for a value of parameter [0;-1;0.5]

res.xopt.feval([0;-1;0.5],'primal')