n-dim state space

[HiroIshida]

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.

[lucattycord]

Hi，thanks for your code！
I want to modify it to an 11-dimensional state space, but there are errors in the number of dimensions. If I want to use the library to write, will I have to modify many parameters?
Thank you so much!

[HiroIshida]

Hi @lucattycord thanks for feedback. I'm currently have trouble in mathworks license, and cannot use MATLAB. I think there needs only slight modification to use in higher dimension case.

What exactly "there are errors in the number of dimensions." is? like, which part of code?

p.s. just few minutes ago, I inquired mathworks about my license stuff, so hopefully I can use matlab in few days and can fix this problem

[lucattycord]

I only modify A, B, C, Q, R, and W matrixes in example_tube MPC.m, and it mentioned me errors above.
The state is 111, and I choose the dim of W_vertex 112. But by calculating, in line 56 of DisturbanceLinearSystem.m, the dim of obj.Ak is 1111 while the dim of obj.W.A is 52 so that it occurs the error.

[HiroIshida]

Could you provide the code (or A, B, C, Q, R values) to reproduce the bug?

Although current compute_mrpi_set is a implementation of this paper, just for ad-hoc remedey, you can use straightforward approximation of invariant set. Cause the definition of the robust invariant set is infinite sum of mikowski addition Z = W ⨁ Ak*W ⨁ Ak^2*W..., you can easily approximate it just by truncation at some point. This approximation, however, leads to Z_approx ⊂ Z, which means Z_approx is not disturbance invariant. So, we must multiply it by some parameter alpha and get Z_approx = alpha*(W ⨁ Ak*W...⨁ Ak^n*W) so that Z ⊂ Z_approx is guaranteed.

However, I cannot find this naive method can work well in higher dimension.

Also, could be nice if you made me a pull request once you can fix the bug.

[lucattycord]

addpath('../src/')
addpath('../src/utils/')

% fix random seed
rng(0);

k_e=10^(-5);
Jx=2562e-6;
Jy=2562e-6;
Jz=5124e-6;
m_uav=1.2954;
r=[0.1;0.1;0.01;3];
input=4;
output=4;
state=11;
T=0.1;%sample time
Np=20;%predictive horizon
Nc=8;%control horizon

% xm(k+1)=A_m*xm(k)+B_m**u(k)

%A_m
A=eye(11,11);
A(1,4)=T;
A(2,5)=T;
A(3,6)=T;
A(4,9)=T/Jx;
A(5,10)=T/Jy;
A(6,11)=T/Jz;
A(7,8)=T+1;
A(9,9)=-10T+1;
A(10,10)=-10T+1;
A(11,11)=-10T+1;
%B_m
B=zeros(11,4);
B(8,1)=T/m_uav;
B(9,2)=10k_eT/Jx;
B(10,3)=10k_eT/Jy;
B(11,4)=10k_e*T/Jz;
%C_m
C=zeros(4,11);
C(1,1)=1;
C(2,2)=1;
C(3,3)=1;
C(4,7)=1;

% make your own discrete linear system with disturbance

Q = diag([1,1,1,1,1,1,1,1,1,1,1]);
R = diag([0.1,0.1,0.1,0.1]);

W_vertex = [0.15,-0.12; 0.1,-0.13; 0.11,-0.11; 0.12,-0.1; 0.14,-0.11; 0.13,-0.14; 0.08,-0.15; 0.09,-0.09; 0.05,0.08; 0.15,-0.1; 0.13,-0.2];
% construct a convex set of disturbance (2dim here)
W = Polyhedron(W_vertex);

% construct disturbance Linear system
disturbance_system = DisturbanceLinearSystem(A, B, Q, R, W);

% constraints on state Xc and input Uc
Xc_vertex = [2, -2; 2 2; -10 2; -10 -2];
Uc_vertex = [1; -1];
Xc = Polyhedron(Xc_vertex);
Uc = Polyhedron(Uc_vertex);

% create a tube_mpc simulater
% if N_horizon is too small, the path will never reach inside the robust MPI-set X_mpi_robust in time step N_horizon, then the problem becomes infeasible.
N_horizon = 10;
w_min = [0; -0.10];
w_max = [0; 0.10];
mpc = TubeModelPredictiveControl(disturbance_system, Xc, Uc, N_horizon);

% The robust MPC guidances the path inside the robust MPI-set so that the path will reach the robust MPI-set in N_horizon.
x = [-7; -2];
savedir_name = 'results';
mkdir(savedir_name);

for i = 1:15
disp(i)
u_next = mpc.solve(x);
x = disturbance_system.propagate(x, u_next); % additive disturbance is considered inside the method
mpc.show_prediction();
saveas(gcf, strcat(savedir_name, '/tmpc_seq', number2string(i), '.png')); % removing this line makes the code much faster
clf;
end

[lucattycord]

My main question lies in the dimension of W. If my status is 11x1 dimension, then due to the upper and lower limits, W is 11x2.

[lucattycord]

Hi,I have modified it to 6-dimensional, you can go to my homepage to see.