H2T Example: Serially connected masses

Alberto Guiggiani - Aug 2013 alberto.guiggiani@imtlucca.it

     >>> Tutorial available online at > http://dysco.imtlucca.it/h2t/masses.html

Physical parameters
System model
LSmodel construction
H2T object + sim data + control parameters
Open-loop simulation
MPT Controller
PnPMPC Controller

This script builds the experiment where M objects of mass m are serially connected by springs and dampeners and attached between two walls. Actuators exert tensions between consecutive masses trying to regulate their positions and velocities to the origin.

MPT3, PnPMPC and WIDE toolbox controllers are tested. System inputs are actuator tensions, while states are masses positions followed by masses velocities. System is state-feedback.

     *** Be sure to run h2t_init before! ***

MPT3, PnPMPC and WIDE toolbox required to run the respective sections.

Run

         >> showdemo Masses_start

to open help, or

         >> echodemo Masses_start

for an interactive demo.

Physical parameters

clear, clc

M = 3;          % Number of masses
m = 1;          % Mass
k = 1.5;         % Spring coefficient
c = 1;          % Damping coefficient

System model

nx = 2*M;       % Number of states
nu = M-1;       % Number of inputs

% Continous-time matrices
Ac = full([zeros(M,M) eye(M)
     gallery('tridiag',M,k/m,-2*k/m,k/m) gallery('tridiag',M,c/m,-c*k/m,c/m)]);

Bc = full([zeros(M,M-1);gallery('tridiag',M-1,-1/m,1/m,0);zeros(1,M-2) -1/m]);

Cc = eye(nx);

sysC = ss(Ac, Bc, Cc, 0);

% Sampling time
Ts = 0.5;

% Discretization
sys = c2d(sysC, Ts);

% Signal names
for i = 1:nu
    sys.InputName(i) = { strcat('actuator', int2str(i)) };
end

for i = 1:M
    sys.StateName(i) = { strcat('pos', int2str(i)) };
    sys.OutputName(i) = { strcat('pos', int2str(i)) };
end

for i = (M+1):nx
    sys.StateName(i) = { strcat('vel', int2str(i-M)) };
    sys.OutputName(i) = { strcat('vel', int2str(i-M)) };
end

LSmodel construction

model = LSmodel(sys, sys.InputName, sys.StateName);

% Signal limits
xl = -4;
xu = 4;
ul = -1;
uu = 1;

for i = 1:nu
    model = model.set_sig_lim(strcat('actuator', int2str(i)), ul, uu);
end

for i = 1:M
    model = model.set_sig_lim(strcat('pos', int2str(i)), xl, xu);
end

for i = M+1:nx
    model = model.set_sig_lim(strcat('vel', int2str(i-M)), xl, xu);
end

model

Large scale model (total order=6, 1 subsystem(s), 0 summator(s)):
  --- Subsystems -----------
    M1: 2 inputs, 6 outputs, order=6, name="M1"
  --- External Inputs -----------
    IN01: actuator1        X  (min=-1,max=1)
    IN02: actuator2        X  (min=-1,max=1)
  --- External Outputs ----------
    OUT01: pos1            X  (min=-4,max=4)
    OUT02: pos2            X  (min=-4,max=4)
    OUT03: pos3            X  (min=-4,max=4)
    OUT04: vel1            X  (min=-4,max=4)
    OUT05: vel2            X  (min=-4,max=4)
    OUT06: vel3            X  (min=-4,max=4)

H2T object + sim data + control parameters

h2t_obj = h2t(model);

% Set prediction horizon
h2t_obj.setParameters('N', 10);

% Set weights
h2t_obj.setParameters('InputWeights', 0.1, 'StateWeights', 1);

h2t_obj

Prediction horizon updated with N=10.000000.

Input weights updated.

State weights updated.

h2t_obj = 

  h2t with properties:

         weights: [1x1 struct]
              nx: 6
              nu: 2
              ny: 6
               N: 10
              Ts: 0.5000
             Net: [2x6 logical]
              X0: []
              Mc: []
               k: []
               m: []
               p: []
     sdpsettings: []
    installedTBX: {'WIDE'  'PnPMPC'  'MPT'}
               M: [1x1 LSmodel]

Open-loop simulation

Tsim = 50;
x = zeros(nx,Tsim);
x(:,1) = [3 -2 1 0 0 0]';
u = zeros(nu,Tsim);

% Open-loop evolution
for i = 2:Tsim
    x(:,i) = sys.A*x(:,i-1) + sys.B*u(:,i-1);
end

% Plot
figure(1), subplot(211)
hold on
plot(x(1,:)), plot(x(2,:),'r'), plot(x(3,:),'k')
title('Open-loop response'), xlabel('Time'), ylabel('Masses displacement')

subplot(212)
hold on
plot(u(1,:)), plot(u(2,:),'r')
title('Control inputs'), xlabel('Time'), ylabel('Exerted force')

% Save sim results
openloop_data.x = x;
openloop_data.u = u;

MPT Controller

% Check parameters for MPT controller
h2t_obj.checkParameters('MPT')

% Build MPT controller
ctrlMPT = h2t_obj.buildController('MPT');

Tsim = 50;
x = zeros(6,Tsim);
x(:,1) = [3 -2 1 0 0 0]';
u = zeros(2,Tsim);

% Evolution with MPT controller
for i = 2:Tsim
    u(:,i-1) = ctrlMPT.evaluate(x(:,i-1));
    x(:,i) = sys.A*x(:,i-1) + sys.B*u(:,i-1);
end

% Plot
figure(2), subplot(211)
hold on
plot(x(1,:)), plot(x(2,:),'r'), plot(x(3,:),'k')
title('MPT Controller'), xlabel('Time'), ylabel('Masses displacement')

subplot(212)
hold on
plot(u(1,:)), plot(u(2,:),'r')
title('Control inputs'), xlabel('Time'), ylabel('Exerted force')

% Save sim results
mpt_data.x = x;
mpt_data.u = u;

------------------------------

Checking options for MPT toolbox...

For details type:
	help MPCController

Required:
	 - Sampling time: OK
	 - Prediction horizon: OK

Facultative:
	 - Weights: OK

------------------------------

------------------------------

Building MPT controller...

SUCCESS

------------------------------

PnPMPC Controller

* Be sure that "extra" folder of PnPMPC toolbox is in Matlab Path *

% Check parameters for PnPMPC controller
h2t_obj.checkParameters('PnPMPC')

% Set missing parameters (for an explanation of their meaning type 'help createCtrlPnPMPC4lsmodel')
h2t_obj.setParameters('k', {10}, 'm', ones(nu,1) ,'p', ones(nx,1), 'sdpsettings', sdpsettings('verbose', 0))

% Build PnPMPC controller
%                               >>> WARNING <<<
% If you get a java exception, simply run the next command again
%                               >>> WARNING <<<
ctrlPnPMPC = h2t_obj.buildController('PnPMPC');
ctrlPnPMPC = ctrlPnPMPC.zeroTerminal(eye(nx), eye(nu));

Tsim = 50;
x = zeros(6,Tsim);
x(:,1) = [3 -2 1 0 0 0]';
u = zeros(2,Tsim);

% Evolution with PnPMPC controller
for i = 2:Tsim
    u(:,i-1) = ctrlPnPMPC.uRH(x(:,i-1),[],[],[],sdpsettings('verbose',0));
    x(:,i) = sys.A*x(:,i-1) + sys.B*u(:,i-1);
end

figure(3), subplot(211)
hold on
plot(x(1,:)), plot(x(2,:),'r'), plot(x(3,:),'k')
title('PnPMPC Controller'), xlabel('Time'), ylabel('Masses displacement')

subplot(212)
hold on
plot(u(1,:)), plot(u(2,:),'r')
title('Control inputs'), xlabel('Time'), ylabel('Exerted force')

% Save sim results
pnpmpc_data.x = x;
pnpmpc_data.u = u;

------------------------------

Checking options for PnPMPC toolbox...

For details type:
	help createCtrlPnPMPC4lsmodel

Required:
	 - Prediction horizon: OK
	 - k parameter: MISSING!
		 Use .setParameters('k',Value)
	 - Input types: MISSING!
		 Use .setParameters('m',Value)
	 - Output types: MISSING!
		 Use .setParameters('p',Value)

Facultative:
	 - YALMIP sdpsettings: MISSING!
		 Use .setParameters('sdpsettings',Value)
		 (type 'help sdpsettings' for infos)
------------------------------

k parameter updated.

Input types updated.

Output types updated.

Sdpsettings updated.

------------------------------

Building PnPMPC controller...
LSS: Time continuous system, discretize with sampletime= 5.000000e-01
CREATECTRLPNP: Controller M1 done

SUCCESS

------------------------------