## Contents

- Set bisection limits for lam.
- Initialize a PDE structure.
- Initialize settings for solving the LPI
- Perform bisection on the value of lam

    % 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.