Introduction

DEMO1_Simple_Stability_Simulation_and_Control_Problem.m See Chapter 2 of the manual for a description.

This document illustrates, with a simple example, how PIETOOLS can be used to simulate systems dynamics, analyse stability and design optimal controllers for the system. The example is the damped wave equation with dynamic boundary input: xi_{tt}=c xi_{ss}-b xi_{t}+sw(t) using the state phi=(X1,X2) such that X1=xi_{s} and X2=xi_{t}: ODE x_{t} = -x + u(t) PDEs: X1_{t} = X2_{s} X2_{t} = cX1_{s}-bX2_{t}+sw(t) With BCs X1(s=0) = 0 X2(s=1) = x = exp(-t)x(t=0)+ int_0^t exp(-t+tau) u(tau) d tau

and regulated output z= int_0^1 xi_{s} ds = xi(s=1)-xi(s=0).

Contents

First, we clear the workspace of any interfering variables

clear all; clc;close all;

Declaring variables for symbolic manipulation

Define the parameters and variables: Here we declare symbolic variables t (st) and s (sx) that stand for time and space, respectively, in polynomial (symbolic) object format. Then, we define states {phi, x}, inputs {w, u}, and output {z}, which will be used to define the PDE system described earlier.

c=1;b=.1;
pvar t s; syms st sx;
phi=state('pde',2);x=state('ode');
w= state('in');u=state('in');
z=state('out',2);

Create the system

Here we use a sys() class object to store the equations that describe the PDE model. First, we initialize the object and then equations are added to the object using addequation() function. Finally, setControl() function is used designate u as a control input.

odepde = sys();
eq_dyn = [diff(x,t,1)==-x+u
                diff(phi,t,1)==[0 1; c 0]*diff(phi,s,1)+[0;s]*w+[0 0;0 -b]*phi];
eq_out= z ==[int([1 0]*phi,s,[0,1])
                                 u];
odepde = addequation(odepde,[eq_dyn;eq_out]);
odepde= setControl(odepde,[u]); % set the control signal
bc1 = [0 1]*subs(phi,s,0)==0; % add the boundary conditions
bc2 = [1 0]*subs(phi,s,1)==x;
odepde = addequation(odepde,[bc1;bc2]);
Initialized sys() object of type "pde"
5 equations were added to sys() object
1 inputs were designated as controlled inputs
2 equations were added to sys() object

Simulating the Open Loop system

Here we simulate the PDE system without any control input to see the open-loop behaviour of the system using PIESIM. The following options are used to define simulation parameters, initial conditions and disturbance inputs for the simulation.

opts.plot = 'no';   % Do not plot the final solution
opts.N = 8;         % Expand using 8 Chebyshev polynomials
opts.tf = 10;        % Simulate up to t = 1;
opts.dt = 1e-2;     % Use time step of 10^-3
opts.intScheme=1;   % Time-step using Backward Differentiation Formula (BDF)
ndiff = [0,2,0];    % The PDE state involves 2 first order differentiable state variables

% To simulate the solution to the PDE system without controller we can use
% the following code.
uinput.ic.PDE = [0,0];
uinput.ic.ODE = 0;
uinput.u=0;
uinput.w = sin(5*st)*exp(-st);
[solution,grids] = PIESIM(odepde, opts, uinput, ndiff);

% Extract actual solution at each time step and defining discretized variables.
tval = solution.timedep.dtime;
phi1 = reshape(solution.timedep.pde(:,1,:),opts.N+1,[]);
phi2 = reshape(solution.timedep.pde(:,2,:),opts.N+1,[]);
zval =solution.timedep.regulated;
wval=subs(uinput.w,st,tval);

% Plots of open-loop system.
figure(1);
surf(tval,grids.phys,phi2,'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row');
h=colorbar ;
colormap jet
box on
ylabel(h,'$|\dot{\mathbf{x}}(t,s)|$','interpreter', 'latex','FontSize',15)
set(gcf, 'Color', 'w');
xlabel('$t$','FontSize',15,'Interpreter','latex');    ylabel('$s$','FontSize',15,'Interpreter','latex');
zlabel('$\dot{\mathbf{x}}(t,s)$','FontSize',15,'Interpreter','latex');
title('Open loop zero-state response with $w(t)=sin(5t)e^{-t}$','Interpreter','latex','FontSize',15);

figure(2);
plot(tval,wval,'k',tval,zval(1,:),'r','LineWidth',2)
grid on
box on
set(gcf, 'Color', 'w');
legend('$\mathbf{w}(t)$','$\mathbf{r}(t)$','Interpreter','latex','FontSize',15)
xlabel('$t$','FontSize',15,'Interpreter','latex');
ylabel('$\mathbf{r}(t)$','FontSize',15,'Interpreter','latex');
title('Open loop zero-state response with $w(t)=sin(5t)e^{-t}$','Interpreter','latex','FontSize',15);
Warning: option  Norder is not defined. Setting to a default value of 2
Solving PDE problem

Encountered 2 state components: 
 x₁(t),      of size 1, finite-dimensional;
 x₂(t,s),    of size 2, differentiable up to order (1) in variables (s);

Encountered 1 actuator input: 
 u(t),    of size 1;

Encountered 1 exogenous input: 
 w(t),    of size 1;

Encountered 1 regulated output: 
 z(t),    of size 2;

Encountered 2 boundary conditions: 
 F₁(t) = 0, of size 1;
 F₂(t) = 0, of size 1;

The order of the state components x has not changed.
Watning: uinput.ifexact is not specified. Defaulting  to false

 --- Reordering the state components to allow for representation as PIE ---

The order of the state components x has not changed.

 --- Converting ODE-PDE to PIE --- 
Setting up Chebyshev matrices for the PIE system
Setup completed: integrating in time
Time integration scheme is numerically stable for the given problem.
Any observed instabilities must be physical.
Solution of an ODE state 1 at a final time 10.000000 is   0.0000
Value of a regulated output 1 at a final time 10.000000 is   0.0146
Value of a regulated output 2 at a final time 10.000000 is   0.0000

Stability Analysis of the system

To use LPIs for the test of internal stability, we have to first convert the PDE representation to a PIE representation. This is performed using the following code.

PIE = convert(odepde,'pie');
PIE = PIE.params;
T = PIE.T;
A = PIE.A;      C1 = PIE.C1;    B2 = PIE.B2;
B1 = PIE.B1;    D11 = PIE.D11;  D12 = PIE.D12;

% Pick parameters for the LPI optimization problem and perform stability
% test.
settings = lpisettings('light');
[prog, P] = lpisolve(PIE,settings,'stability');

% Other LPI tests:
% Similarly, we can search for an input-to-output L2-gain bound using LPIs
% as shown below.
[prog, P, gamma] = lpisolve(PIE,settings,'l2gain');
 --- Reordering the state components to allow for representation as PIE ---

The order of the state components x has not changed.

 --- Converting ODE-PDE to PIE --- 
Initialized sys() object of type "pde"
Conversion to pie was successful

 --- Executing Primal Stability Test --- 
- Parameterizing Positive Lyapunov Operator using specified options...
- Constructing the Negativity Constraint...
- Enforcing the Negativity Constraint...
  - Using an Equality constraint...
- Solving the LPI using the specified SDP solver...
Size: 755  181
 
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 181, order n = 44, dim = 756, blocks = 4
nnz(A) = 4607 + 0, nnz(ADA) = 32585, nnz(L) = 16383
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            9.60E-01 0.000
  1 :  -1.99E-06 2.66E-01 0.000 0.2773 0.9000 0.9000   1.00  1  1  1.5E+00
  2 :   4.50E-06 9.03E-02 0.000 0.3393 0.9000 0.9000   1.00  1  1  5.2E-01
  3 :   1.62E-06 3.27E-02 0.000 0.3621 0.9000 0.9000   1.00  1  1  1.9E-01
  4 :   1.79E-07 8.69E-03 0.000 0.2657 0.9000 0.9000   1.00  1  1  5.0E-02
  5 :   8.34E-10 2.20E-04 0.000 0.0253 0.9900 0.9900   1.00  1  1  1.3E-03
  6 :  -4.24E-14 2.44E-09 0.000 0.0000 1.0000 1.0000   1.00  1  1  1.8E-08
  7 :  -1.81E-17 2.07E-12 0.000 0.0008 0.9999 0.9996   1.00  1  1  1.2E-11

iter seconds digits       c*x               b*y
  7      0.4  13.9  0.0000000000e+00 -1.8097152689e-17
|Ax-b| =   1.8e-11, [Ay-c]_+ =   1.5E-12, |x|=  1.2e+01, |y|=  6.6e+00

Detailed timing (sec)
   Pre          IPM          Post
1.880E-01    4.060E-01    3.100E-02    
Max-norms: ||b||=2.000000e-04, ||c|| = 0,
Cholesky |add|=0, |skip| = 94, ||L.L|| = 115.464.
 
Residual norm: 1.8292e-11
 
         iter: 7
    feasratio: 1.0000
         pinf: 0
         dinf: 0
       numerr: 0
       timing: [0.1880 0.4060 0.0310]
      wallsec: 0.6250
       cpusec: 0.5625


 --- Searching for Hinf gain bound using primal KYP lemma --- 
- Declaring Positive Lyapunov Operator variable using specified options...
- Constructing the Negativity Constraint...
- Enforcing the Negativity Constraint...
  - Using an Equality constraint...
- Solving the LPI using the specified SDP solver...
Size: 931  265
 
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Detected 1 diagonal SDP block(s) with 1 linear variables
Put 1 free variables in a quadratic cone
eqs m = 265, order n = 53, dim = 933, blocks = 5
nnz(A) = 5490 + 0, nnz(ADA) = 66439, nnz(L) = 33352
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            5.03E-01 0.000
  1 :   1.45E+00 1.32E-01 0.000 0.2618 0.9000 0.9000   0.60  1  1  1.9E+00
  2 :   1.99E+00 4.73E-02 0.000 0.3589 0.9000 0.9000   1.01  1  1  6.9E-01
  3 :   3.00E+00 1.49E-02 0.000 0.3162 0.9000 0.9000   0.55  1  1  3.0E-01
  4 :   3.90E+00 4.72E-03 0.000 0.3158 0.9000 0.9000   0.61  1  1  1.1E-01
  5 :   4.52E+00 1.90E-03 0.000 0.4027 0.9000 0.9000   0.86  1  1  4.8E-02
  6 :   4.94E+00 6.29E-04 0.000 0.3310 0.9000 0.9000   1.12  1  1  1.5E-02
  7 :   5.11E+00 1.73E-04 0.000 0.2748 0.9000 0.9000   1.14  1  1  3.8E-03
  8 :   5.15E+00 4.83E-05 0.000 0.2797 0.9000 0.9000   1.07  1  1  1.0E-03
  9 :   5.16E+00 1.55E-05 0.000 0.3211 0.9000 0.9000   1.02  1  1  3.3E-04
 10 :   5.17E+00 1.75E-06 0.000 0.1128 0.9000 0.8977   0.94  1  1  1.0E-04
 11 :   5.17E+00 5.62E-07 0.000 0.3210 0.8686 0.9000   0.85  2  2  3.7E-05
 12 :   5.17E+00 2.34E-07 0.000 0.4167 0.9000 0.9109   0.76  2  2  1.7E-05
 13 :   5.17E+00 1.04E-07 0.000 0.4440 0.9000 0.9139   0.69  2  2  8.6E-06
 14 :   5.17E+00 4.68E-08 0.000 0.4504 0.9000 0.9216   0.64  5  5  4.5E-06
 15 :   5.17E+00 2.38E-08 0.000 0.5083 0.9000 0.9384   0.61  9  6  2.6E-06
 16 :   5.17E+00 1.43E-08 0.000 0.6027 0.9000 0.8950   0.65  9  9  1.8E-06
 17 :   5.17E+00 8.55E-09 0.254 0.5962 0.9000 0.9342   0.65 12  9  1.2E-06
 18 :   5.17E+00 5.50E-09 0.311 0.6437 0.9000 0.9316   0.64 20 19  8.1E-07
Run into numerical problems.

iter seconds digits       c*x               b*y
 18      0.4   6.0  5.1664118557e+00  5.1664172384e+00
|Ax-b| =   4.0e-06, [Ay-c]_+ =   4.0E-08, |x|=  2.4e+02, |y|=  1.5e+04

Detailed timing (sec)
   Pre          IPM          Post
3.100E-02    3.930E-01    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=6, |skip| = 131, ||L.L|| = 5.24206e+06.
 
Residual norm: 3.9914e-06
 
         iter: 18
    feasratio: 0.6360
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0.0310 0.3930 0]
      wallsec: 0.4240
       cpusec: 0.4688

The H-infty norm of the given system is upper bounded by:
    5.1664

Find Hinf-optimal controller

Next, we wish to find a controller that improves the input-to-output performance. For that, we can solve the Hinf-optimal controller LPI using the following code. The function, if the optimization problem is successfully solved, returns the controller, Kval and the L2-gain metric for the closed loop system, gam_val.

[prog, Kval, gam_val] = lpisolve(PIE, settings,'hinf-controller');
 --- Executing Search for H_infty Optimal Controller --- 
- Declaring Positive Storage Operator variable and indefinite Controller operator variable using specified options...
- Parameterize the derivative inequality...
- Enforcing the Negativity Constraint...
  - Using an Equality constraint...
- Solving the LPI using the specified SDP solver...
Size: 938  263
 
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Detected 1 diagonal SDP block(s) with 1 linear variables
Put 8 free variables in a quadratic cone
eqs m = 263, order n = 53, dim = 940, blocks = 5
nnz(A) = 4601 + 0, nnz(ADA) = 65061, nnz(L) = 32662
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            6.20E-01 0.000
  1 :   8.10E-01 1.66E-01 0.000 0.2683 0.9000 0.9000   1.20  1  1  1.9E+00
  2 :   9.58E-01 6.12E-02 0.000 0.3680 0.9000 0.9000   1.55  1  1  5.7E-01
  3 :   9.56E-01 2.01E-02 0.000 0.3280 0.9000 0.9000   1.24  1  1  1.7E-01
  4 :   8.80E-01 6.56E-03 0.000 0.3268 0.9000 0.9000   1.12  1  1  5.4E-02
  5 :   8.30E-01 2.43E-03 0.000 0.3701 0.9000 0.9000   1.08  1  1  2.0E-02
  6 :   8.05E-01 8.86E-04 0.000 0.3651 0.9000 0.9000   0.97  1  1  7.4E-03
  7 :   8.05E-01 1.57E-04 0.152 0.1767 0.9000 0.0000   0.87  1  1  4.1E-03
  8 :   7.95E-01 2.71E-05 0.000 0.1728 0.9270 0.9000   0.88  1  1  9.5E-04
  9 :   7.86E-01 7.94E-06 0.000 0.2934 0.9067 0.9000   0.69  1  1  3.2E-04
 10 :   7.83E-01 2.87E-06 0.000 0.3615 0.6270 0.9000   0.77  1  1  1.4E-04
 11 :   7.81E-01 1.30E-06 0.000 0.4544 0.9000 0.6939   0.75  1  1  7.1E-05
 12 :   7.80E-01 4.64E-07 0.000 0.3557 0.9000 0.9000   0.82  1  1  2.7E-05
 13 :   7.79E-01 1.49E-07 0.000 0.3222 0.9000 0.9000   0.84  1  1  9.5E-06
 14 :   7.79E-01 4.08E-08 0.000 0.2733 0.9000 0.9000   0.89  1  1  2.8E-06
 15 :   7.79E-01 1.12E-08 0.000 0.2736 0.9000 0.9000   0.91  1  1  7.9E-07
 16 :   7.79E-01 2.34E-09 0.000 0.2096 0.9000 0.9000   0.95  1  1  1.7E-07
 17 :   7.79E-01 1.72E-10 0.112 0.0735 0.9900 0.9900   0.98  1  2  1.3E-08
 18 :   7.79E-01 1.39E-11 0.201 0.0809 0.9900 0.9900   1.00  2  2  1.0E-09

iter seconds digits       c*x               b*y
 18      0.3   8.7  7.7865323963e-01  7.7865324107e-01
|Ax-b| =   5.8e-10, [Ay-c]_+ =   5.9E-10, |x|=  4.5e+01, |y|=  3.4e+00

Detailed timing (sec)
   Pre          IPM          Post
1.600E-02    2.340E-01    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=0, |skip| = 132, ||L.L|| = 3897.34.
 
Residual norm: 5.7884e-10
 
         iter: 18
    feasratio: 0.9977
         pinf: 0
         dinf: 0
       numerr: 0
       timing: [0.0160 0.2340 0]
      wallsec: 0.2500
       cpusec: 0.2969

The closed-loop H-infty norm of the given system is upper bounded by:
    0.7787

Constructing closed-loop system and simulation

Now, we construct the closed-loop system using the controller obtained by solving the above LPI and then re-runs the simulations to see if there was any improvement in the performance.

PIE_CL = closedLoopPIE(PIE,Kval);
PIE_CL = pie_struct(PIE_CL);
PIE_CL = initialize(PIE_CL);

% Simulate the solution to the PIE with controller
[solution_CL,grids] = PIESIM(PIE_CL,opts,uinput,ndiff);
Warning: option  Norder is not defined. Setting to a default value of 2
Solving PIE problem
Setting up Chebyshev matrices for the PIE system
Setup completed: integrating in time
Time integration scheme is numerically stable for the given problem.
Any observed instabilities must be physical.
Solution of an ODE state 1 at a final time 10.000000 is   0.0027
Value of a regulated output 1 at a final time 10.000000 is  -0.0020
Value of a regulated output 2 at a final time 10.000000 is  -0.0011

Plotting the closed loop response

Having found an optimal controller, constructed closed loop system, and running PIESIM to simulate, we can now extract the solution and plot it against the open-loop response to see the benefit in using the controller. Firstly, we notice that the L2-gain bound has significantly lowered (from 5.2 to 0.78). Then, looking at the output in presence of a disturbance, we see that the effect of disturbance on this neutrally stable system is reduced.

tval = solution_CL.timedep.dtime;
phi1 = reshape(solution_CL.timedep.pde(:,1,:),opts.N+1,[]);
phi2 = reshape(solution_CL.timedep.pde(:,2,:),opts.N+1,[]);
zval_cl =solution_CL.timedep.regulated;
wval=subs(uinput.w,st,tval);

% Plots Closed Loop.
figure(3)
surf(tval,grids.phys,phi2,'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row');
h=colorbar ;
colormap jet
box on
ylabel(h,'$|\dot{\mathbf{x}}(t,s)|$','interpreter', 'latex','FontSize',15)
set(gcf, 'Color', 'w');
xlabel('$t$','FontSize',15,'Interpreter','latex');    ylabel('$s$','FontSize',15,'Interpreter','latex');
zlabel('$\dot{\mathbf{x}}(t,s)$','FontSize',15,'Interpreter','latex');
title('Closed loop zero-state response with $w(t)=sin(5t)e^{-t}$','Interpreter','latex','FontSize',15);

figure(4);
plot(tval,wval,'k',tval,zval_cl(1,:),'r',tval,zval_cl(2,:),'b','LineWidth',2)
grid on
box on
set(gcf, 'Color', 'w');
legend('$\mathbf{w}(t)$','$\mathbf{r}(t)$','$\mathbf{u}(t)$','Interpreter','latex','FontSize',15)
xlabel('$t$','FontSize',15,'Interpreter','latex');
ylabel('$\mathbf{r}(t)$','FontSize',15,'Interpreter','latex');
title('Closed loop zero-state response with $w(t)=sin(5t)e^{-t}$','Interpreter','latex','FontSize',15);

figure(5);
plot(tval,zval(1,:),'r--',tval,zval_cl(1,:),'r','LineWidth',2)
grid on
box on
set(gcf, 'Color', 'w');
legend('$\mathbf{r}(t)$','$\mathbf{r}_{cl}(t)$','Interpreter','latex','FontSize',15)
xlabel('$t$','FontSize',15,'Interpreter','latex');
ylabel('$\mathbf{z}(t)$','FontSize',15,'Interpreter','latex');
title('Zero-state response with $w(t)=sin(5t)e^{-t}$','Interpreter','latex','FontSize',15);