Contents

% DEMO4_stability_parameter_bisection.m
% See Chapter 11.4 of the manual for a description.
%
% This document illustrates how PIETOOLS can be used to find the maximal
% value of a parameter for which a PDE is stable, using bisection.

% Reaction-diffusion PDE
% \dot{x}(t,s) = lam*x(t,s) + x_{ss}(t,s);
% x(t,0) = x(t,1) = 0;
% Stable when lam <= pi^2 = 9.8696 (Ahmadi 2015).
clc; clear;
echo on
%%%%%%%%%%%%%%%%%% Start Code Snippet %%%%%%%%%%%%%%%%%%

Set bisection limits for lam.

lam_min = 0;        lam_max = 20;
lam = 0.5*(lam_min + lam_max);
n_iters = 8;

Initialize a PDE structure.

a = 0;  b = 1;
pvar s
pde_struct PDE;
PDE.x{1}.vars = s;
PDE.x{1}.dom = [a,b];

% Set the PDE \dot{x}(t,s) = lam*x(t,s) + x_{ss}(t,s);
PDE.x{1}.term{1}.C = lam;
PDE.x{1}.term{2}.D = 2;

% Set the BCs x(t,a) = x(t,b) = 0;
PDE.BC{1}.term{1}.loc = a;
PDE.BC{2}.term{1}.loc = b;

Initialize settings for solving the LPI

settings = lpisettings('veryheavy');
if strcmp(settings.sos_opts.solver,'sedumi')
    settings.sos_opts.params.fid = 0;   % Suppress output in command window
end

Perform bisection on the value of lam

for iter = 1:n_iters
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;

    % Update the PIE.
    PIE = convert(PDE,'pie');

    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])
    prog = PIETOOLS_stability(PIE,settings);

    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
        % Stability cannot be verified, decreasing the value of lam...
        lam_max = lam;
        lam = 0.5*(lam_min + lam_max);
    else
        % The system is stable, trying a larger value of lam...
        lam_min = lam;
        lam = 0.5*(lam_min + lam_max);
    end
end
%%%%%%%%%%%%%%%%%% End Code Snippet %%%%%%%%%%%%%%%%%%
echo off

fprintf(['\n Stability of the system could be verified for lam<=',num2str(lam_min),'.\n'])
fprintf([' An analytic bound on lam guaranteeing stability is pi^2=',num2str(pi^2),'.\n']);

% @article{valmorbida2015stability,
%   title={Stability analysis for a class of partial differential equations via semidefinite programming},
%   author={Valmorbida, Giorgio and Ahmadi, Mohamadreza and Papachristodoulou, Antonis},
%   journal={IEEE Transactions on Automatic Control},
%   volume={61},
%   number={6},
%   pages={1649--1654},
%   year={2015},
%   publisher={IEEE}
% }
%%%%%%%%%%%%%%%%%% Start Code Snippet %%%%%%%%%%%%%%%%%%

%%% Set bisection limits for lam.
lam_min = 0;        lam_max = 20;
lam = 0.5*(lam_min + lam_max);
n_iters = 8;

%%% Initialize a PDE structure.
a = 0;  b = 1;
pvar s
pde_struct PDE;
PDE.x{1}.vars = s;
PDE.x{1}.dom = [a,b];

% Set the PDE \dot{x}(t,s) = lam*x(t,s) + x_{ss}(t,s);
PDE.x{1}.term{1}.C = lam;
PDE.x{1}.term{2}.D = 2;

% Set the BCs x(t,a) = x(t,b) = 0;
PDE.BC{1}.term{1}.loc = a;
PDE.BC{2}.term{1}.loc = b;

%%% Initialize settings for solving the LPI
settings = lpisettings('veryheavy');
if strcmp(settings.sos_opts.solver,'sedumi')
    settings.sos_opts.params.fid = 0;   % Suppress output in command window
end

%%% Perform bisection on the value of lam
for iter = 1:n_iters
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 10 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 1.4202e-06
 
         iter: 17
    feasratio: -0.2636
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0.0310 0.5290 0.0110]
      wallsec: 0.5710
       cpusec: 1.0156


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
        % Stability cannot be verified, decreasing the value of lam...
        lam_max = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 5 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 1.2274e-07
 
         iter: 13
    feasratio: 1.0001
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0 0.3600 0]
      wallsec: 0.3600
       cpusec: 0.6562


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
    else
        % The system is stable, trying a larger value of lam...
        lam_min = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 7.5 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 1.8976e-06
 
         iter: 12
    feasratio: 1.0041
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0 0.3050 0]
      wallsec: 0.3050
       cpusec: 0.4688


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
    else
        % The system is stable, trying a larger value of lam...
        lam_min = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 8.75 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 2.8761e-07
 
         iter: 14
    feasratio: 0.9876
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0.0070 0.3350 0]
      wallsec: 0.3420
       cpusec: 0.3906


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
    else
        % The system is stable, trying a larger value of lam...
        lam_min = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 9.375 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 2.1765e-07
 
         iter: 15
    feasratio: 1.0077
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0 0.4200 0.0160]
      wallsec: 0.4360
       cpusec: 0.6875


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
    else
        % The system is stable, trying a larger value of lam...
        lam_min = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 9.6875 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 4.3523e-07
 
         iter: 15
    feasratio: 1.0099
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0.0070 0.3440 0]
      wallsec: 0.3510
       cpusec: 0.4375


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
    else
        % The system is stable, trying a larger value of lam...
        lam_min = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 9.8438 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 3.0045e-07
 
         iter: 16
    feasratio: 1.0122
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0 0.4790 0.0010]
      wallsec: 0.4800
       cpusec: 0.8750


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
    else
        % The system is stable, trying a larger value of lam...
        lam_min = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
    % Update the value of lam in the PDE.
    PDE.x{1}.term{1}.C = lam;
    
    % Update the PIE.
    PIE = convert(PDE,'pie');

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

 --- Converting ODE-PDE to PIE --- 
    
    % Re-run the stability test.
    fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])

 --- Running the stability test with lam = 9.9219 ---
    prog = PIETOOLS_stability(PIE,settings);

 --- 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: 3122   230
 
 
Residual norm: 7.3406e-07
 
         iter: 17
    feasratio: 0.1010
         pinf: 0
         dinf: 0
       numerr: 1
       timing: [0 0.3920 0]
      wallsec: 0.3920
       cpusec: 0.4375


    % Check if the system is stable
    if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
        % Stability cannot be verified, decreasing the value of lam...
        lam_max = lam;
        lam = 0.5*(lam_min + lam_max);
    end    
end
%%%%%%%%%%%%%%%%%% End Code Snippet %%%%%%%%%%%%%%%%%%
echo off

 Stability of the system could be verified for lam<=9.8438.
 An analytic bound on lam guaranteeing stability is pi^2=9.8696.