LowCom - holaza/mptplus GitHub Wiki
Design and verification of low-complexity explicit MPC controllers in MPT3
Reference:
The features of LowCom package were presented in:
M. Kvasnica – J. Holaza – B. Takács – D. Ingole: Design and Verification of Low-Complexity Explicit MPC Controllers in MPT3. In Proceedings of the European Control Conference 2015, Linz, Austria, pp. 2600–2605, doi.org/10.1109/ECC.2015.7330929, 2015.
The polynomial approximation of the designed Explicit MPC piece-wise affine (PWA) control law is implemented according to: M. Kvasnica – J. Löfberg – M. Fikar: Stabilizing polynomial approximation of explicit MPC. Automatica, vol. 10 (47), pp. 2292–2297, 2011.
Low-Complexity Explicit MPC
Let us consider the following model of the system, along with the corresponding constraints, penalty settings, and the prediction horizon:
model = LTISystem('A', [1 1; 0 1], 'B', [1; 0.5])
model.x.min = [-5; -5], model.x.max = [5, 5]
model.u.min = -1, model.u.max = 1
model.x.penalty = QuadFunction(eye(2))
model.u.penalty = QuadFunction(1)
model.x.with('terminalPenalty')
model.x.terminalPenalty = QuadFunction(5*eye(2))
model.x.with('terminalSet')
model.x.terminalSet = Polyhedron('lb',[-1;-1], 'ub', [1;1])
N = 5
The explicit model predictive controller can be obtained when running:
eMPC = MPCController(model,N).toExplicit()
Now we are able to plot the polytopic partition using the following command:
eMPC.partition.plot()
Complexity Reduction without Inducing Suboptimality
First, we introduce the set of complexity reduction methods without inducing the suboptimality into the evaluation of the control law.
Optimal-region-merging-based complexity reduction
The first complexity reduction method is based on merging together critical regions which share the same expression of the feedback law, and whose union is convex. The method is described in more detail in:
The optimal-region-merging (ORM) technique is achieved by using the simplify function as follows:
eMPCsimple = eMPC.simplify('orm')
which returns a new simplified explicit controller object which behaves in the same way as the original object eMPC. The ORM method is applicable to arbitrary explicit MPC feedback laws regardless of their properties (e.g. discontinuity of the feedback law). The corresponding simplified polytopic partition can be plotted analogously to the original explicit MPC:
eMPCsimple.partition.plot()
Clipping-based complexity reduction
The clipping-based method exploits the continuity of the explicit MPC feedback law to remove critical regions where the control action is saturated at minimal or maximal value. The space covered by such saturated regions is then replaced by extensions of the unsaturated critical regions, followed by applying a clipping filter. The technique can be found in:
This complexity reduction can be achieved by:
eMPCsimple = eMPC.simplify('clipping')
Separation-based complexity reduction
Lastly, the separation-based technique exploits saturation properties of continuous explicit MPC feedback laws. In contrast to the clipping method, the saturated regions are removed and only the unsaturated critical regions are kept. Then a separation function is devised to uniquely determine the value of the optimal control input for initial conditions which are not contained in the retained part of the explicit solution. The complexity reduction is described in
and can be invoked using the following command:
eMPCsimple = eMPC.simplify('separation')
Complexity Reduction with Suboptimality
Next, we introduce the set of complexity reduction methods by inducing the suboptimality into the evaluation of the control law.
Fitting-based complexity reduction
The two-step fitting procedure was suggested in:
In the first step, a simpler partition composed of fewer critical regions is designed by solving the MPC problem for a shorter value of the prediction horizon. Subsequently, the parameters of the simpler feedback are optimized to minimize the suboptimality represented by the integrated squared error between the original and the simpler feedback. Such a complexity reduction scheme is available via:
eMPCsimple = eMPC.simplify('fitting')
The average and worst-case decrease of performance of the simplified controller versus the optimal one can be compared by:
[mean_performance, worst_case_performance] = eMPC.comparePerformance(eMPCsimple)
Minimum-time-based complexity reduction
The second method to reduce complexity by introducing suboptimality is suggested in:
In this strategy, the complexity is reduced by solving a simpler cost function. Here, a minimum-time problem is solved where the objective is to minimize the arrival time to a certain terminal set. The controller can be simplified by:
eMPCsimple = EMinTimeController(model)
Polynomial-approximation-based complexity reduction
The polynomial approximation constructed for the piece-wise affine (PWA) control law determined by solving the Explicit MPC with a linear (1-, Inf- norm) cost function according to: M. Kvasnica – J. Löfberg – M. Fikar: Stabilizing polynomial approximation of explicit MPC. Automatica, vol. 10 (47), pp. 2292–2297, 2011.
The approximated polynomial controller can be simplified, evaluated, and plotted by typing:
eMPCsimple = PolynomialMPC(EMPC)
u = eMPCsimple.evaluate(x0)
eMPCsimple.plot
Demo:
%% Polynomial Approximation
% LTI system
model = LTISystem('A', [1.2], 'B', [1]);
model.u.min = [ -1 ];
model.u.max = [ 1 ];
model.x.min = [ -4 ];
model.x.max = [ 4 ];
% Linear penalty functions
model.x.penalty = OneNormFunction(diag([ 1 ]));
model.u.penalty = OneNormFunction(diag([ 1 ]));
% Prediction horizon
N = 10;
% Construct MPC controller
iMPC = MPCController(model,N)
% Explicit MPC controller construction
eMPC = iMPC.toExplicit
% MPC evaluation
x0 = [ -2 ]; % Initial condition
[ u_explicit, feasible, openloop ] = eMPC.evaluate(x0) % Explicit MPC evaluation
% Polynomial Approximation Settings
% degree of the approximation polynomial
polyMPC = PolynomialMPC(eMPC)
polyMPC.plot
[ u_approx, feasible ] = polyMPC.evaluate(x0)
Behind the Scenes
In this section, we discuss the basic building blocks of the aforementioned techniques. The basic functionality required by the ORM method is the recognition of the convexity of a union of polyhedra. Given polytope P1 and a polytope P2. To determine, whether their union is convex, one can use the following method:
P1 = Polyhedron(H1, h1)
P2 = Polyhedron(H2, h2)
convex = isConvex(PolyUnion([P1, P2]))
More specifically, the answer is obtained by checking whether the difference set of the convex hull of the union and the union itself is an empty set. The convex hull Ph can be obtained by enumerating the vertices V of the polytopes P1 and P2 and by removing the redundant ones. The isConvex method, therefore, corresponds to the following procedure:
V = [P1.V; P2.V] % vertex enumeration
Ph = Polyhedron(V) % convex hull
diff = Ph \ [P1, P2] % set difference
convex = diff.isEmptySet() % comparison
The clipping and saturation-based complexity reduction techniques rely on the identification of critical regions of the PWA feedback law in which the control action is saturated. The structure containing the information about the saturated and unsaturated regions can be obtained by:
eMPC.optimizer.findSaturated('primal')
Another tool provides the ability to construct a linear separator p of nonconvex sets S1 and S2 given as unions of polytopes Pi a and Qi. The separator's parameters can be found by:
S1 = PolyUnion([P1, ..., Pn])
S2 = PolyUnion([Q1, ..., Qm])
p = SeparationController.separate(S1,S2)
In the fitting-based approach, the main technical difficulty is represented by the integration of a function f over polytope P in generic n-dimensional Euclidean space. Such a procedure can be achieved with the following steps:
P = Polyhedron('lb', [-1; -2], 'ub', [3; 4])
f = QuadFunction([1 0.1; 0.1 2], [1 0], -3.5)
P.addFunction(f, 'fx')
int = P.integrate('fx')
Another functionality that can be beneficial is the tessellation of polytope P into simplices. It can be achieved by:
D = P.triangulate()
The volume of every i-th simplex can then be calculated via:
vol = D(i).volume()
Verification and analysis
The purpose of the verification and analysis is to verify whether a closed-loop system, composed of a discrete-time system and an MPC controller, exhibits certain safety properties. Let us consider the LTI system described with model2. The closed-loop dynamics is then created by:
model2 = LTIsystem('A', [1 1; 0 0.9], 'B', [1; 0.5])
loop = ClosedLoop(eMPC, model2)
Such a closed-loop system can then be analyzed e.g. by performing simulations with the initial state x0 and Nsim simulation steps:
data = loop.simulate(x0, Nsim)
A more rigorous way to verify closed-loop stability for all feasible initial conditions is to devise a Lyapunov function. We are able to synthesize piecewise quadratic (PWQ) and piecewise affine (PWA) Lyapunov functions if the system’s dynamics and the explicit representation of the MPC feedback law are available. In such a case the closed-loop system becomes a PWA system and can be obtained from:
pwa = loop.toSystem()
The Lyapunov function for the autonomous PWA system is then computed by:
lyap = pwa.lyapunov(type)
The type can be specified subject to the required form of the function, i.e., two options (pwq for piece-wise-quadratic, pwa for piece-wise-affine, ) are available:
type='pwq'
type='pwa'
Such an approach can be also used to verify closed-loop stability when a suboptimal controller simple is employed. We are able to verify that both the fitting controller as well as the minimum-time strategy described previously provide closed-loop stability by running:
pwa = ClosedLoop(eMPCsimple, model).tomodeltem()
lyap = pwa.lyapunov('pwq')
It is also possible to verify the invariance of a given closed-loop system. The invariance of an autonomous PWA system can be assessed by:
answer = pwa.isInvariant()
If the answer is negative, we can calculate the invariant subset via:
pwa_inv = pwa.invariantSet()
The computation of invariant sets is in fact more general since it covers the construction of maximal control invariant sets as well. If the system’s dynamics is linear of piecewise affine, the maximal control invariant set is constructed by:
Cinf = model.invariantSet()
In the calculation procedure of the maximal control invariant sets, the backward reachability analysis is performed. This can be achieved based on the recursion algorithm, based on the initial set S(k), while the pre-set S(k+1) can be calculated by:
S(k+1) = model.reachableSet('X', S(k), 'direction', 'backward')
One case where maximal invariant sets are used frequently is when the LQR set, where constraints are not active, needs to be constructed. Specifically, the LQR set is obtained via:
OLQR = model.LQRSet()
It is also possible to compute the forward reachable sets starting from some initial set S via:
F = model.reachableSet('X', S, 'direction', 'forward')
Depending on the level of granularity the user wants to operate with, basis reachability tasks can be performed by applying direct and inverse affine transformations of polyhedra. If A and c would represent matrices of an autonomous affine system x(k+1) = Ax(k) + c, then the direct transform yields forward reachable sets, while the inverse transform coincides with backward reachability. The transformations are available via the * operator:
direct = A*S+c
inverse = S*A+c