Contents
clc; clear;
echo on
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];
PDE.x{1}.term{1}.C = lam;
PDE.x{1}.term{2}.D = 2;
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;
end
Perform bisection on the value of lam
for iter = 1:n_iters
PDE.x{1}.term{1}.C = lam;
PIE = convert(PDE,'pie');
fprintf(['\n',' --- Running the stability test with lam = ',num2str(lam),' ---\n'])
prog = PIETOOLS_stability(PIE,settings);
if prog.solinfo.info.dinf || prog.solinfo.info.pinf || abs(prog.solinfo.info.feasratio - 1)>0.3
lam_max = lam;
lam = 0.5*(lam_min + lam_max);
else
lam_min = lam;
lam = 0.5*(lam_min + lam_max);
end
end
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']);
%%%%%%%%%%%%%%%%%% 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.