hiroishida/robust-tube-mpc

An example code for robust model predictive control using tube

Robust Model Predictive Control Using Tube

This repository includes examples for the tube model predictive control (tube-MPC)[1] as well as the generic model predictive control (MPC) written in MATLAB.

Requirement

1. optimization_toolbox (matlab)
   
2. control_toolbox (matlab)
   
3. Multi-Parametric Toolbox 3 (open-source and freely available at http://people.ee.ethz.ch/~mpt/3/)

Feedback, bug reports, contributions

If you find this package helpful, giving a "star" to this repositry will be a happy feedback for me! If you find a bug, or have more broader kind of quession about tube MPC,please post that in the issue page. I will try hard to respond to questions via e-mail but, I strongly recommend do it in the issue page. It's much easier for me to keep myself on track.

Usage

See example/example_tubeMPC.m and example/example_MPC.m for the tube-MPC and generic MPC, respectively. Note that every inequality constraint here is expressed as a convex set. For example, the constraints on the state Xc is specified as a rectangular, which is constructed with 4 vertexes. When considering a 1-dim input Uc, Uc will be specified by min and max value (i.e. u∊[u_min, u_max]), so it will be constructed by 2 vertexes. For more detail, please see the example codes.

Short introduction to the tube MPC

After running example/example_tubeMPC.m, you will get the following figure sequence.

Now that you can see that the green nominal trajectory starting from the bottom left of the figure and surrounding a "tube". At each time step, the nominal trajectory (green line) is computed online.

Let me give some important details. The red region Xc that contains the pink region Xc-Z is the state constraint that we give first. However, considering the uncertainty, the tube-MPC designs the nominal trajectory to be located inside Xc-Z, which enables to put "tube" around the nominal trajectory such that the tube is also contained in Xc-Z. Of course, the input sequence associated with the nominal trajectory is inside of Uc-KZ.

Disturbance invariant set

I think one may get stuck at computation of what paper [1] called "disturbance invariant set". The disturbance invariant set is an infinite Minkowski addition Z = W ⨁ Ak*W ⨁ Ak^2*W..., where ⨁ denotes Minkowski addition. Because it's an infinite sum of Minkowski addition, computing Z analytically is intractable. In [2], Racovic proposed a method to efficiently compute an outer approximiation of Z, which seems to be heavily used in MPC community. In this repository, computation of Z takes place in the constructor of DisturbanceLinearSystem class. To understand how Z guarantee the robustness, running example/example_dist_inv_set.m may help you.

Maximum positively invariant set

I used the maximal positively invariant (MPI) set Xmpi as the terminal constraint set. (Terminal constraint is usually denoted as Xf in literature). Book [3] explains the concept of the MPI and algorithm well in section 2.4. Xmpi is computed in the constructor of OptimalControler.m. Note that the MPI set is computed with Xc and Uc in the normal MPC setting, but in the tube-MPC the MPI set is computed with Xc⊖Zand Uc⊖Z instead.

Reference

[1] Mayne, David Q., María M. Seron, and S. V. Raković. "Robust model predictive control of constrained linear systems with bounded disturbances." Automatica 41.2 (2005): 219-224. [2] Rakovic, Sasa V., et al. "Invariant approximations of the minimal robust positively invariant set." IEEE Transactions on Automatic Control 50.3 (2005): 406-410. [3] Kouvaritakis, Basil, and Mark Cannon. "Model predictive control." Switzerland: Springer International Publishing (2016).

robust-tube-mpc's Issues

How to track a time-varying reference signal?

How to track a time-varying reference signal?

n-dim state space

Current version support only 2-dim state space. It will be better if the software can deal with arbitrary dim state space. To this end, there will be two tasks required:
Task1) support n-dim space in OptimalControler.m.
Task2) develop a visualization tool which can handle n-dim state space. [if possible]

Task2 might be difficult as some nice idea is necessary to neatly visualize over 4dim spaces.

the caculation of mRPI set

hello! I change your model to three order model，
where A = [1，0.2，-1; 0， 1， -0.2;0，0， 0.6],
B= [ 0; 0; 0.6],
Q = diag([1, 1,0.1]),
R = 0.1 ;
W_vertex = [5 5 5;5 5 -5;5 -5 5;5 -5 -5;-5 5 5;-5 5 -5;-5 -5 5;-5 -5 -5];
During the caculation of mRPI set of disturbance linear system, this setting leads to the value of alpha converges to epsilon/(epsilon + Ms) very slowly ,so the value of s is very big. Then the calculation of Fs is beyond the ability of MATLAB. Could you please give some guidance? Thank you very much.

Code requirements

An error "Error using mskoptimset" occurs inside quadprog

Error using mskoptimset (line 67)
Expected argument 2 to be a string parameter name or an options structure
created with mskoptimset.

Error in quadprog (line 78)
options = mskoptimset(defaultopt,options);

Error in OptimalControler/solve (line 76)
[var_optim, ~, exitflag] = quadprog(obj.H, [], C_ineq1_relaxed, C_ineq2, C_eq1, C_eq2, [], [], [], options);

Error in ModelPredictiveControl/solve (line 21)
[x_nominal_seq, u_nominal_seq] = obj.optcon.solve();

Error in example_MPC (line 27)
u_next = mpc.solve(x);

Do you know how I can solve this problem? Thank you very much!

Robust output feedback model predictive control of constrained linear systems仿真

你好，请问有没有这篇论文的仿真Robust output feedback model predictive control ofconstrained linear systems，小白刚开始有点不懂，求指教

How to extend to distributed mpc case

Hi, thank you for your code again!

Can this code extend to distributed model predictive control? I see that the OptimalControler class is only for objective function V(s) = s'HS. While in the distributed form, the cost function is V(s) = sum(x_i-x_j)'Q(x_i-x_j)+u_i'Ru_i for all its neighbors j and only x_i is the decision variable, x_j is known. So, do you know how to update this code for the dsitributed cost function?

Thank you very much!

question about initial constraints

Thank you for sharing the code. I have a doubt about constraints on initial variable in problem P*(x). In equation (21) of reference [1] (p.222), the initial constraints are expressed as "x belongs to x0+Z ", where x0 is the decision variable of optimal control problem P*(x), and x is the known current state. While in the file "TubeModelPredictiveControl.m", line 41-42, the constraints are described as "Xinit = x_init + obj.sys.Z", where "x_init" is the known current state x, and "Xinit" is the constraint boundary added to the problem with the form Ax<=b to restrict initial decision variable. I think that the code in line 41 implies "x0 belongs to x+Z", which is not consistent with the original constraint "x belongs to x0+Z ". What's your opinion?

Error in control input for larger noise levels

Hi!

First of all, thanks for this repository! I think it is very useful for anyone who is starting to get in touch with tube-based MPC.

I was testing out the example in the repository and notice something peculiar. When I set the noise level from 0.15 to 0.2, i.e.,
from:
W_vertex = [0.15, 0.15; 0.15, -0.15; -0.15, -0.15; -0.15, 0.15];
to:
vertex_max = 0.2;
W_vertex = [vertex_max, vertex_max; vertex_max, -vertex_max; -vertex_max, -vertex_max; -vertex_max, vertex_max];

I observe that
(i) the values of the final control input (u0 + K(x-x0)) exceed the constraint set (Uc) in the first timestep, and
(ii) the values of the nominal control input (u0) exceed the robust constraint set (Uc_robust) in the first timestep.

Here is a plot of the time histories of the final and nominal control inputs. The blue line is the final control input and the red line is the nominal control input. The control limits are at [-1, 1].

I think we should expect the final control input to stay within Uc, is that correct?

Hope to hear from you and happy to discuss. Thanks!

Sophisticated way to generate disturbance invariant set

Currently, disturbance invariant set Z is generated by truncation of infinite Minkowski addition: Z = W ⨁ Ak*W ⨁ Ak^2*W.. as mentioned in README.

It would be better if this Z can be generated by the more sophisticated way such as the one of Rakovic [1].

[1] Rakovic, Sasa V., et al. "Invariant approximations of the minimal robust positively invariant set." IEEE Transactions on Automatic Control 50.3 (2005): 406-410.

caculation of mRPI set

Hi,

I have some doubts about the calculation of mRPI set. In equation (11) of ref[2], there is a transpose of $A^s$, however, in line 56, there isn't a transpose of obj.Ak^s. Is there any reason for the inconsistency?

Similar, in equation(13) of ref[2], there is a transpose of $A^i$, but there isn't in line 59.

Thanks

Error using quadprog

From a user, I got a report a error similar to https://groups.google.com/forum/#!topic/mpt-user/SL7bYNsKNTQ. This error may occur in running example_tubeMPC.

According to the user the following way will be a workaround:
Using ‘interior-point-convex’ instead of ’trust-region-reflective’ seems working now. I just want to report this problem and hope other people having the same problem will benefit from my experience. The Matlab version from my side is 2019a.

The author (me) currently couldn't reproduce this error, but I will investigate it later. I tested only on ubuntu16 + matlab2020a.

Refactoring is needed

As I wrote this code when I had just started to study the object oriented programming, the code is messy especially in the following point: The current code uses inheritance in the wrong way. All of inheritances (is-a) in the code must be replaced by compositions (has-a).

quadprog usage

Thanks for sharing your nice code.

Could you please explain why you use the quadprog command in a loop in OptimalController.m) and what is the role of the following command? (in optimalcontroler.m line 79):

C_ineq1_relaxed = C_ineq1_relaxed*0.999;

Best Regard

Trying to understand how F and G matrices are constructed

I am trying to understand, how F and G matrices are constructed from X and U matrixes

How have you decided these are the constraints for U and X: https://github.com/HiroIshida/robust-tube-mpc/blob/master/example/example_tubeMPC.m#L19

Thanks

some confusion

hi, appreciating your code on RMPC,
is there possible to revise it into a model with 4 state+1 input+2output ?

The "tube" can be feasible ? hoping for answer

Error

Hello,
I had a few questions
Please guide me.

The first is that it gives this error when I run your example code.

Undefined function or variable 'Polyhedron'.

Error in example_dist_inv_set (line 12)
W = Polyhedron(W_vertex);

Best regards,

questions: running examples

Hi,

I have been trying to run the examples. However, I have not been able to do it. The main issue is about the command "mpc.simulate(Tsimu, x_init)" in the examples. I have got the following errors:
Error using mskoptimset (line 67)
Expected argument 2 to be a string parameter name or an options structure
created with mskoptimset.
Error in quadprog (line 78)
options = mskoptimset(defaultopt,options);
Error in OptimalControler/solve (line 62)
[var_optim, ~, exitflag] = quadprog(obj.H, [], C_neq1_relaxed, obj.C_neq2, obj.C_eq1, obj.C_eq2(x_init), [], [], [], options);
Error in TubeModelPredictiveControl/simulate (line 52)
[x_nominal_seq, u_nominal_seq] = obj.optcon.solve(x_init);

Do you know why I have got these errors?

result of example.m (disturbance invariant MPC)

question on computing mrpi

Hi,

I have a 4-order linear discreted system with A, B, Q, R, when I use the function compute_mrpi_set() to solve the maximum robust invariant set, the result of Miniski addition is not convergered. I mean the alpha is becoming larger with the increasing iterative steps, while the correct alpha should become smaller. I want to know what's the reason, it is caused by the feedback maricx Ak or disturbance polyhedron W？

Your answer would be appreciated.

Yongxiang