% Example used in the paper submitted at ROCOND 2015
% Title: R-RoMulOC: A unified tool for randomized and robust multiobjective control
% Authors: Mohammadreza Chamanbaz, Fabrizio Dabbene, Dimitri Peaucelle, Roberto Tempo  

clear all
uncertainty('clear')
%% Fisrt part for defining the models

% Mechanical system described by:
%    d^2x     dx 
%  M ---- + D -- + K x = E w
%    dt^2     dt
%
%  z = C x + F w
%
% The K and F matrices are known:
n=3;m=3;p=3;
K = [  3  -1  0 ;
  0   7  1 ;
  -1   0  4 ];
F = [ 0 0 0 ;
  0 0 0 ;
  1 0 0 ];
% The  inertia M is uncertain such that M = M0 + M1*Dm*M2 with
M0 = [  1  1  0 ;
  -1  1  1 ;
  0 -1  1 ];
M1 = [ 1 ;
  0 ;
  -1 ];
M2 = [ 1 0  0 ;
  0 0 -1 ];
% and Dm is a 1-by-2 time-varying matrix with values evolving in some
% ellispoïd
eli = melli([-1 0.5 ; 0.5 -1],[2;-2],5/0.25);
% Next we define the uncertainty as in the
% ellipsoid, the last argument is just a Tag
Dm = udiss(eli, 'Uncert on Inertia');
% The damping D matrix is uncertain such that D = D0 + D1*Dd*D2 with
D0 = [ 0.5 0   0   ;
  0   0.9 0   ;
  0   0   0.7 ];
D1 = [ 1 0 ;
  1 1 ;
  1 0 ];
D2 = [ 1 0 0 ;
  0 0 1 ];
% and Dd is a 2-by-2 norm bounded parametric matrix: Dd''*Dd<=0.25
Dd = unb( 2, 2, 0.25, 'Uncert on Damping');
% The input matrix E is uncertain such that E = E0 + E1*De*E2 with\n');
E0 = eye(3);
E1 = [ 1 ; 0 ; 0 ];
E2 = [ 0 0 1 ];
% and De is a scalar in the interval [-0.25 0.25]
De = uinter(-1, 1, 'Uncertt on Input');
% The output matrix C is uncertain such that C = C0 + C1*Dc*C2 with
C0 = eye(3);
C1 = [ 0 1 ;
  1 0 ;
  0 0 ];
C2 = [ 1 0 0 ;
  0 1 0 ];
% and Dc is a 2-by-2 matrix in the convex hull of the three matrices (in
% that special case of polytopes with three vertices it is possible to draw
% samples uniformly in the set)
Dcv(:,:,1)=[-0.25,0;0,-0.25];
Dcv(:,:,2)=[ 0.25,0.25;0,-0.25];
Dcv(:,:,3)=[ 0.25,0.25;0, 0.25];
Dc = upoly( Dcv, 'Uncertainties on the Output');
% This uncertain system can be modeled as an LFT-type uncertain system.
% That is: the feedback connection of the uncertain operator
% such that     wd = Delta*zd    with 
Delta = diag(Dm,Dd,De,Dc);
% and an LTI system such that
%
% dX/dt = A  X + Bd  wd + Bw  w
%  zd   = Cd X + Dd  wd + Ddw w
%  z    = Cz X + Dzd wd + Dzw w
%
% We shall not explain how these matrices are obtained but it is easy to see
% by hand that it does correspond to the problem
sys=ssmodel('Demo Example 1');
iM0=inv(M0);
sys.A = [ zeros(n) , eye(n) ; -iM0*D0  , -iM0*K ];
sys.Bd((n+1):2*n,:)=[ -iM0*M1 , -iM0*D1 , iM0*E1 , zeros(n,size(C1,2)) ];
sys.Bw((n+1):2*n,:) = -iM0*E0 ;
sys.Cd = [-M2*iM0*D0 , -M2*iM0*K ;...
    zeros(size(D2,1),n) , D2 ;...
    zeros(size(E2,1),2*n) ;...
    C2 , zeros(size(C2,1),n) ];
sys.Ddd(1:size(M2,1),:) = [ -M2*iM0*M1 , -M2*iM0*D1 , M2*iM0*E1 , zeros(size(M2,1),size(C1,2)) ];
sys.Ddw = [ M2*iM0*E0 ; zeros(size(D2,1),m) ; E2; zeros(size(C2,1),m) ];
sys.Cz = [C0 , zeros(p,n) ];
sys.Dzd(:,(end-size(C1,2)+1):end) = C1;
sys.Dzw = F;
% As in the Romuloc Demo, the input/output signals used to define the
% performance are shaped 
g=tf([0.1 1],[0.02 1]);
G(1,1)=g;G(2,2)=g;G(3,2)=0;
G=ssmodel(G,'filter');
sys = shape( sys, G);
sys = shape( [1 0 0], sys, 'shaped model');
% we add some input signals for the state-feedback problem
sys.Bu=sys.Bw;
sys.Ddu=sys.Ddw;
% All the data of the uncertain system is finally gathered in a USSMODEL object:
usys=ussmodel(sys,Delta)
% When looking at the display of this system you see it is an LFT
% (inevitable because there is a need to model the inverse of the uncertain
% inertia of the system). 

%% How one can draw samples
usample(usys.uncertainty)
usample(usys)

%% Randomized performance computation
quiz_Pr = ctrpb('analysis', 'rand')+hinfty(usys); 
opts=randsettings('epsilon',1e-3,'delta',1e-6); 
H_Pr=solvesdp(quiz_Pr,opts)

%% Randomized Performance Verification
quiz_Pr2 = ctrpb('analysis', 'rand')+hinfty(usys, 2.498);
opts=randsettings('epsilon',1e-2,'delta',1e-2);
solvesdp(quiz_Pr2,opts)

%% Deterministic Worst-Case Performance
quiz_pilf = ctrpb('analysis', 'PILF')+hinfty(usys);
H_pilf=solvesdp(quiz_pilf)

%% Improved Deterministic Worst-Case Performance
quiz_pdlf = ctrpb('analysis', 'PDLF')+hinfty(usys);
H_pdlf=solvesdp(quiz_pdlf)

%% Deterministic multiobjective design
quiz_sf_d = ctrpb('state-feedback')+hinfty(usys,1)+dstability(usys,region('disc',0,10));
K_d=solvesdp(quiz_sf_d)
%% Randomized scenario approach
quiz_sf_s = ctrpb('state-feedback','rand')...
    +hinfty(usys,1)+dstability(usys,region('disc',0,10));
opts=randsettings('epsilon',1e-1,'delta',1e-2,'method','scenario');
K_s=solvesdp(quiz_sf_s,opts)

%% Randomized sequential approaches
quiz_sf_g = ctrpb('state-feedback','rand')...;
    +hinfty(usys,1)+dstability(usys,region('disc',0,10));
opts=randsettings('epsilon',1e-1,'delta',1e-2,'method','gradient','Nout',1e5);
K_g=solvesdp(quiz_sf_g,opts)

