## Contents

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%

    % DEMO5_Hinf_optimal_estimator.m
    % See Chapter 11.5 of the manual for a description.
    %
    % This document illustrates how an Hinfty optimal estimator can be designed
    % for a PIE, and how the estimated state can be simulated.
    %
    % We consider the system defined by:
    % PDE:      \dot{x}(t,s) = (d^2/ds^2) x(t,s) + lam x(t,s) + w(t),   s in [0,1];
    % Outputs:       z(t)    = int_{0}^{1} x(t,s) ds + w(t);
    %                y(t)    = x(t,1);
    % BCs:                0  = x(t,0) = x_{s}(t,1);
    % (unstable for lam > pi^2/4 = 2.4674)
    %
    % Letting v:=(d^2/ds^2)x, We derive an equivalent PIE of the form:
    % [T \dot{v}](t,s) = [A v](t,s) + [B1 w](t,s);
    %             z(t) = [C1 v](t)  + [D11 w](t);
    %             y(t) = [C2 v](t)  + [D21 w](t);
    %
    % So that x = T*v. We design an estimator of the form:
    % [T \dot{vhat}](t,s) = [A vhat](t,s) + [L(yhat-y)](t,s);
    %             zhat(t) = [C1 vhat](t)
    %             yhat(t) = [C2 vhat](t)
    %
    % Using this estimator, the errors e=(vhat-v) and ztilde = (zhat-z) satisfy
    % [T \dot{e}](t,s) = [(A+L*C2) e](t,s) - [(B1+L*D21) w](t,s);
    %        ztilde(t) = [C1 e](t)         - [D11 w](t);
    %
    % We wish to compute an operator L that minimizes the L2 gain from
    % disturbances w to error ztilde in the output. This is achieved solving
    % the LPI
    %
    % min_{gam,P,Z} gam
    % s.t.  P>=0
    %       [-gam*I,           -D11',  -(P*B1+Z*D21)'*T           ]=: Q <=0
    %       [-D11,             -gam*I, C1                         ]
    %       [-T'*(P*B1+Z*D21), C1',    (P*A+Z*C2)'*T+T'*(P*A+Z*C2)]
    %
    % Then, using L = P^{-1}*Z, the L2 gain satisfied
    % ||ztilde||_{L2}/||w||_{L2} <= gam
    %
    % We manually declare this LPI here, but it can also be solved using the
    % "PIETOOLS_Hinf_estimator" executive file.
    % We simulate the PDE state x and estimated state xhat using PIESIM.
    %

    clc; clear; close all;
    echo on
    %%%%%%%%%%%%%%%%%% Start Code Snippet %%%%%%%%%%%%%%%%%%

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%

% % Declare the PDE, and convert it to a PIE. Declare the PDE using command line parser

    pvar s t                                        % Initialize polynomial variables t (time) and s (space)
    lam = 4;                                        % Set reaction parameter
    PDE = sys();                                    % Initialize the PDE structure
    x = state('pde');   w = state('in');            % Initialize PDE state x and disturbance w
    y = state('out');   z = state('out');           % Initialize outputs y and z
    eqs = [diff(x,t) == diff(x,s,2) + lam*x + w;    % define the PDE equation
           z == int(x,s,[0,1]) + w;                 % define the regulated output equation
           y == subs(x,s,1);                        % define the observed output equation
           subs(x,s,0) == 0;                        % define the first boundary condition
           subs(diff(x,s),s,1) == 0];               % define the second boundary condition
    PDE = addequation(PDE,eqs);                     % Add the equations to the PDE structure
    PDE = setObserve(PDE,y);                        % Set y as an observed output
    display_PDE(PDE);                               % Display the system in the command window

    % Compute the associated PIE.
    PIE = convert(PDE,'pie');       PIE = PIE.params;
    % Extract the PI operators defining the PIE.
    T = PIE.T;
    A = PIE.A;      C1 = PIE.C1;    C2 = PIE.C2;
    B1 = PIE.B1;    D11 = PIE.D11;  D21 = PIE.D21;
    % A PIE with state v, disturbance w, regulated output z, and observed
    % output y, with no disturbance in the BCs, has the structure
    % T * dot{v}(t) = A*v(t)  + B1*w(t)
    %         z(t)  = C1*v(t) + D11*w(t)
    %         y(t)  = C2*v(t) + D12*w(t)

    %%%%%%%%%%%%%%%%%% Start Code Snippet %%%%%%%%%%%%%%%%%%
    %% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
    % % % Declare the PDE, and convert it to a PIE.
    % Declare the PDE using command line parser
    pvar s t                                        % Initialize polynomial variables t (time) and s (space)
    lam = 4;                                        % Set reaction parameter
    PDE = sys();                                    % Initialize the PDE structure
    Initialized sys() object of type "pde"
    x = state('pde');   w = state('in');            % Initialize PDE state x and disturbance w
    y = state('out');   z = state('out');           % Initialize outputs y and z
    eqs = [diff(x,t) == diff(x,s,2) + lam*x + w;    % define the PDE equation
           z == int(x,s,[0,1]) + w;                 % define the regulated output equation
           y == subs(x,s,1);                        % define the observed output equation
           subs(x,s,0) == 0;                        % define the first boundary condition
           subs(diff(x,s),s,1) == 0];               % define the second boundary condition
    PDE = addequation(PDE,eqs);                     % Add the equations to the PDE structure
    5 equations were added to sys() object
    PDE = setObserve(PDE,y);                        % Set y as an observed output
    1 outputs were designated as observed outputs
    display_PDE(PDE);                               % Display the system in the command window

      ∂ₜ x(t,s) = ∂²ₛ x(t,s) + 4 * x(t,s) + w(t);

      y(t) = x(t,1);

      z(t) = ₀∫¹[x(t,s)]ds + w(t);

      0 = - x(t,0);
      0 = - ∂ₛ x(t,1);

    % Compute the associated PIE.
    PIE = convert(PDE,'pie');       PIE = PIE.params;

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

     --- Converting ODE-PDE to PIE ---
    Initialized sys() object of type "pde"
    Conversion to pie was successful
    % Extract the PI operators defining the PIE.
    T = PIE.T;
    A = PIE.A;      C1 = PIE.C1;    C2 = PIE.C2;
    B1 = PIE.B1;    D11 = PIE.D11;  D21 = PIE.D21;

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%

% % Compute an optimal observer operator L for the PIE.

    use_executive = false;  % <-- set to true to use predefined executive
    if use_executive
        % % Use the predefined Hinf estimator executive function.
        % Declare settings to use: choose from
        % extreme < stripped < light < heavy < veryheavy   or   custom
        % Heavier settings may increase accuracy, but also computation time.
        settings = lpisettings('heavy');
        % Call the Hinf estimator executive:
        % Returns a solved LPI program structure, an optimal observer operator
        % Lval, and a bound gam_val on the L2-gain ||zhat-z||/||w|| for the
        % associated estimator.
        [prog, Lval, gam_val] = PIETOOLS_Hinf_estimator(PIE, settings);
    else
        % % Manually construct and solve the LPI program for optimal
        % % estimator synthesis.

        % Initialize the LPI program
        vars = PIE.vars(:);         % Extract the spatial variables
        prog = sosprogram(vars);    % Initialize an LPI program structure in the considered spatial variables

        % Declare the decision variable gamma
        dpvar gam;
        prog = sosdecvar(prog, gam);

        % Declare a positive semidefinite PI operator decision variable P>=0
        Pdim = T.dim(:,1);              % Row dimensions of the operator P
        Pdom = PIE.dom;                 % Spatial domain of the operator P
        Pdeg = {6,[2,3,5],[2,3,5]};     % Degrees of monomials used to define P (call "help poslpivar" for more info)
        opts.sep = 1;                   % Set P.R.R1=P.R.R2 to reduce computational complexity
        [prog,P] = poslpivar(prog,Pdim,Pdom,Pdeg,opts);
        %eppos = 1e-6;
        %P.R.R0 = P.R.R0 + eppos*eye(size(P));

        % Declare the indefinite PI operator decision variable Z
        Zdim = C2.dim(:,[2,1]);         % Row and column dimensions of the operator Z
        Zdom = PIE.dom;                 % Spatial domain of the operator Z
        Zdeg = [4,0,0];                 % Degrees of monomials defining Z (call "help lpivar" for more info)
        [prog,Z] = lpivar(prog,Zdim,Zdom,Zdeg);

        % Declare the LPI constraint Q<=0.
        nw = size(B1,2);    nz = size(C1,1);
        Q = [-gam*eye(nw),     -D11',        -(P*B1+Z*D21)'*T;
             -D11,             -gam*eye(nz), C1;
             -T'*(P*B1+Z*D21), C1',          (P*A+Z*C2)'*T+T'*(P*A+Z*C2)];
        prog = lpi_ineq(prog,-Q);       % Add the constraint -Q>=0 to the LPI program

        % Set the objective function: minimize gam
        prog = sossetobj(prog, gam);

        % Solve the optimization program
        opts.solver = 'sedumi';         % Use SeDuMi to solve the SDP
        opts.simplify = true;           % Simplify the SDP before solving
        prog_sol = sossolve(prog,opts);

        % Extract the solved value of gam and the operators P and Z
        gam_val = sosgetsol(prog_sol,gam);
        Pval = getsol_lpivar(prog_sol,P);
        Zval = getsol_lpivar(prog_sol,Z);

        % Build the optimal observer operator L.
        Lval = getObserver(Pval,Zval);
    end

    % A PIE with state v, disturbance w, regulated output z, and observed
    % output y, with no disturbance in the BCs, has the structure
    % T * dot{v}(t) = A*v(t)  + B1*w(t)
    %         z(t)  = C1*v(t) + D11*w(t)
    %         y(t)  = C2*v(t) + D12*w(t)

    %% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
    % % % Compute an optimal observer operator L for the PIE.

    use_executive = false;  % <-- set to true to use predefined executive
    if use_executive
    else
        % % Manually construct and solve the LPI program for optimal
        % % estimator synthesis.

        % Initialize the LPI program
        vars = PIE.vars(:);         % Extract the spatial variables
        prog = sosprogram(vars);    % Initialize an LPI program structure in the considered spatial variables

        % Declare the decision variable gamma
        dpvar gam;
        prog = sosdecvar(prog, gam);

        % Declare a positive semidefinite PI operator decision variable P>=0
        Pdim = T.dim(:,1);              % Row dimensions of the operator P
        Pdom = PIE.dom;                 % Spatial domain of the operator P
        Pdeg = {6,[2,3,5],[2,3,5]};     % Degrees of monomials used to define P (call "help poslpivar" for more info)
        opts.sep = 1;                   % Set P.R.R1=P.R.R2 to reduce computational complexity
        [prog,P] = poslpivar(prog,Pdim,Pdom,Pdeg,opts);
        %eppos = 1e-6;
        %P.R.R0 = P.R.R0 + eppos*eye(size(P));

        % Declare the indefinite PI operator decision variable Z
        Zdim = C2.dim(:,[2,1]);         % Row and column dimensions of the operator Z
        Zdom = PIE.dom;                 % Spatial domain of the operator Z
        Zdeg = [4,0,0];                 % Degrees of monomials defining Z (call "help lpivar" for more info)
        [prog,Z] = lpivar(prog,Zdim,Zdom,Zdeg);

        % Declare the LPI constraint Q<=0.
        nw = size(B1,2);    nz = size(C1,1);
        Q = [-gam*eye(nw),     -D11',        -(P*B1+Z*D21)'*T;
             -D11,             -gam*eye(nz), C1;
             -T'*(P*B1+Z*D21), C1',          (P*A+Z*C2)'*T+T'*(P*A+Z*C2)];
        prog = lpi_ineq(prog,-Q);       % Add the constraint -Q>=0 to the LPI program

        % Set the objective function: minimize gam
        prog = sossetobj(prog, gam);

        % Solve the optimization program
        opts.solver = 'sedumi';         % Use SeDuMi to solve the SDP
        opts.simplify = true;           % Simplify the SDP before solving
        prog_sol = sossolve(prog,opts);
    Running simplification process:
    Old A size: 20248  673
    New A size: 19967  672
    Size: 19967    672

    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
    Put 6 free variables in a quadratic cone
    eqs m = 672, order n = 162, dim = 19969, blocks = 4
    nnz(A) = 52030 + 0, nnz(ADA) = 451584, nnz(L) = 226128
     it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
      0 :            8.45E+00 0.000
      1 :   1.58E+00 2.21E+00 0.000 0.2614 0.9000 0.9000   0.51  1  1  2.3E+01
      2 :   1.33E+00 7.15E-01 0.000 0.3234 0.9000 0.9000   1.11  1  1  6.7E+00
      3 :   1.24E+00 2.05E-01 0.000 0.2872 0.9000 0.9000   1.12  1  1  1.9E+00
      4 :   1.17E+00 5.99E-02 0.000 0.2917 0.9000 0.9000   1.05  1  1  5.3E-01
      5 :   1.11E+00 1.91E-02 0.000 0.3197 0.9000 0.9000   0.86  1  1  2.0E-01
      6 :   1.08E+00 7.23E-03 0.000 0.3779 0.9000 0.9000   0.78  1  1  8.9E-02
      7 :   1.06E+00 2.72E-03 0.000 0.3767 0.9000 0.9000   0.81  1  1  3.9E-02
      8 :   1.05E+00 1.26E-03 0.000 0.4636 0.9000 0.9000   0.64  1  1  2.4E-02
      9 :   1.04E+00 5.31E-04 0.000 0.4208 0.9000 0.9000   0.77  1  1  1.2E-02
     10 :   1.03E+00 2.49E-04 0.000 0.4687 0.9000 0.9000   0.45  1  1  8.0E-03
     11 :   1.02E+00 9.93E-05 0.000 0.3985 0.9000 0.9000   0.71  1  1  3.7E-03
     12 :   1.02E+00 4.35E-05 0.000 0.4382 0.9000 0.9000   0.46  1  1  2.3E-03
     13 :   1.01E+00 1.73E-05 0.000 0.3987 0.9000 0.9000   0.68  1  2  1.1E-03
     14 :   1.01E+00 7.81E-06 0.000 0.4502 0.9000 0.9000   0.47  1  2  7.1E-04
     15 :   1.01E+00 3.11E-06 0.000 0.3989 0.9000 0.9000   0.69  2  2  3.3E-04
     16 :   1.01E+00 9.79E-07 0.000 0.3144 0.9000 0.0000   0.45  2  3  2.2E-04
     17 :   1.01E+00 2.89E-07 0.000 0.2957 0.9098 0.9000   0.35  2  2  8.2E-05
     18 :   1.01E+00 8.63E-08 0.000 0.2983 0.9000 0.0000   0.10  4  4  5.4E-05
     19 :   1.00E+00 2.38E-08 0.000 0.2758 0.9095 0.9000   0.02  4  5  1.9E-05
    Run into numerical problems.

    iter seconds digits       c*x               b*y
     19      5.0   3.1  1.0035251762e+00  1.0043628463e+00
    |Ax-b| =   2.7e-05, [Ay-c]_+ =   1.6E-06, |x|=  1.3e+03, |y|=  7.0e+03

    Detailed timing (sec)
       Pre          IPM          Post
    4.700E-02    4.413E+00    0.000E+00
    Max-norms: ||b||=1, ||c|| = 1,
    Cholesky |add|=0, |skip| = 335, ||L.L|| = 1.5284e+07.

    Residual norm: 2.6731e-05

             iter: 19
        feasratio: 0.0232
             pinf: 0
             dinf: 0
           numerr: 1
           timing: [0.0470 4.4130 0]
          wallsec: 4.4600
           cpusec: 5.0938

        % Extract the solved value of gam and the operators P and Z
        gam_val = sosgetsol(prog_sol,gam);
        Pval = getsol_lpivar(prog_sol,P);
        Zval = getsol_lpivar(prog_sol,Z);

        % Build the optimal observer operator L.
        Lval = getObserver(Pval,Zval);
    end

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%

% % Build a PIE modeling the actual state v, and estimated state vhat: [T, 0] [\dot{v}(t) ] = [A, 0 ] [v(t) ] + [B1 ] w(t) [0, T] [\dot{vhat}(t)] [-L*C2, A+L*C2] [vhat(t)] + [L*D21]

                [z(t)   ] = [C1, 0 ] [v   ]           + [D11] w(t)
                [zhat(t)]   [0,  C1] [vhat]             [0  ]

This PIE has the form T_CL * \dot{V}(t) = A_CL * V(t) + B_CL * w(t) Z(t) = C_CL * V(t) + D_CL * w(t)

where V = [v; vhat] and Z = [z; zhat].

    % Construct the operators defining the PIE.
    T_CL = [T, 0*T; 0*T, T];        % use 0*T to define zero-operaotr of same dimensions as T
    A_CL = [A, 0*A; -Lval*C2, A+Lval*C2];   B_CL = [B1; Lval*D21];
    C_CL = [C1, 0*C1; 0*C1, C1];            D_CL = [D11; 0*D11];

    % Declare the PIE.
    PIE_CL = pie_struct();                      % Initialize the PIE structure
    PIE_CL.vars = PIE.vars;                     % Set the spatial variables of the PIE
    PIE_CL.dom = PIE.dom;                       % Set the domain of the spatial variables
    PIE_CL.T = T_CL;                            % Declare the operator T
    PIE_CL.A = A_CL;        PIE_CL.B1 = B_CL;   % Declare the operators A and B1
    PIE_CL.C1 = C_CL;       PIE_CL.D11 = D_CL;  % Declare the operators C1 and D11
    PIE_CL = initialize(PIE_CL);                % Fill in all the blanks in the PIE structure

    %% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
    % % % Build a PIE modeling the actual state v, and estimated state vhat:
    % [T, 0] [\dot{v}(t)   ] = [A,     0     ] [v(t)   ] + [B1   ] w(t)
    % [0, T] [\dot{vhat}(t)]   [-L*C2, A+L*C2] [vhat(t)] + [L*D21]
    %
    %              [z(t)   ] = [C1, 0 ] [v   ]           + [D11] w(t)
    %              [zhat(t)]   [0,  C1] [vhat]             [0  ]
    %
    % This PIE has the form
    % T_CL * \dot{V}(t) = A_CL * V(t) + B_CL * w(t)
    %             Z(t)  = C_CL * V(t) + D_CL * w(t)
    %
    % where V = [v; vhat] and Z = [z; zhat].

    % Construct the operators defining the PIE.
    T_CL = [T, 0*T; 0*T, T];        % use 0*T to define zero-operaotr of same dimensions as T
    A_CL = [A, 0*A; -Lval*C2, A+Lval*C2];   B_CL = [B1; Lval*D21];
    C_CL = [C1, 0*C1; 0*C1, C1];            D_CL = [D11; 0*D11];

    % Declare the PIE.
    PIE_CL = pie_struct();                      % Initialize the PIE structure
    PIE_CL.vars = PIE.vars;                     % Set the spatial variables of the PIE
    PIE_CL.dom = PIE.dom;                       % Set the domain of the spatial variables
    PIE_CL.T = T_CL;                            % Declare the operator T
    PIE_CL.A = A_CL;        PIE_CL.B1 = B_CL;   % Declare the operators A and B1
    PIE_CL.C1 = C_CL;       PIE_CL.D11 = D_CL;  % Declare the operators C1 and D11
    PIE_CL = initialize(PIE_CL);                % Fill in all the blanks in the PIE structure

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%

% % Simulate and plot the actual and estimated PDE state using PIESIM

    % Declare initial conditions for the state components of the PIE
    syms st sx real             % Declare symbolix variables st (time) and sx (space)
    uinput.ic.PDE = [-10*sx;    % Set the actual initial PIE state value
                     0];        % Set the estimated initial PIE state value

    % Declare the value of the disturbance w(t)
    uinput.w = 2*sin(pi*st);

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

    % Simulate the solution to the PIE with estimator.
    [solution,grid] = PIESIM(PIE_CL,opts,uinput,ndiff);

    % Plot the solution at several grid points using the subroutine at the
    % bottom of this DEMO.
    grid_idcs = [1,5,7];        % Only plot at first, fifth and seventh grid points
    plot_figure_Hinf_optimal_estimator_DEMO(solution,grid,opts,grid_idcs)

    %% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
    % % % Simulate and plot the actual and estimated PDE state using PIESIM

    % Declare initial conditions for the state components of the PIE
    syms st sx real             % Declare symbolix variables st (time) and sx (space)
    uinput.ic.PDE = [-10*sx;    % Set the actual initial PIE state value
                     0];        % Set the estimated initial PIE state value

    % Declare the value of the disturbance w(t)
    uinput.w = 2*sin(pi*st);

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

    % Simulate the solution to the PIE with estimator.
    [solution,grid] = PIESIM(PIE_CL,opts,uinput,ndiff);
    Warning: option  Norder is not defined. Setting to a default value of 2
    Solving PIE problem
    Warning: control inputs are not defined. Defaulting to zero
    Initial conditions on ODE states are not defined. Defaulting to 1
    Too many initial conditions are defined. Ignoring extra conditions
    Setting up Chebyshev matrices for the PIE system
    Setup completed: integrating in time
    Time integration scheme is numerically unstable for the given problem.
    Try increasing time step to 2.610
    or decreasing an order of the scheme (opts.Norder).
    Value of a regulated output 1 at a final time 1.000000 is  12.0320
    Value of a regulated output 2 at a final time 1.000000 is  12.0322

    % Plot the solution at several grid points using the subroutine at the
    % bottom of this DEMO.
    grid_idcs = [1,5,7];        % Only plot at first, fifth and seventh grid points
    plot_figure_Hinf_optimal_estimator_DEMO(solution,grid,opts,grid_idcs)

    %%%%%%%%%%%%%%%%%% End Code Snippet %%%%%%%%%%%%%%%%%%
    echo off

    %%
    %%%%%%%%%%%%%%%%%% End Code Snippet %%%%%%%%%%%%%%%%%%
    echo off

## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%

    function plot_figure_Hinf_optimal_estimator_DEMO(solution,grid,opts,grid_idcs)
    % % % Plot simulated values of actual state and estimated state at several
    % % % grid points, as specified by grid_idcs.

    % Extract actual and estimated solution at each time step.
    x_act = reshape(solution.timedep.pde(:,1,:),opts.N+1,[]);
    x_est = reshape(solution.timedep.pde(:,2,:),opts.N+1,[]);
    tval = solution.timedep.dtime;

    % Set options for the plot
    plot_indcs = floor(logspace(0,log(opts.tf/opts.dt)/log(10),40));
    %plot_indcs = floor(linspace(1,opts.tf/opts.dt,66));
    tplot = tval(plot_indcs);           % Only plot at select times
    colors = {'b','g','m','r','k','c','r','y','o'};     % Colors for the plot

    % Plot evolution of actual and estimated
    fig1 = figure(1);
    for j=grid_idcs
        s_pos = num2str(grid.phys(j));  % Position associated to grid index.
        subplot(1,2,1)
        hold on
        plot(tplot,x_act(j,plot_indcs),[colors{j},'-'],'LineWidth',2,'DisplayName',['$\mathbf{x}(s=',s_pos,')$']);
        plot(tplot,x_est(j,plot_indcs),[colors{j},'o--'],'LineWidth',1.5,'DisplayName',['$\mathbf{\hat{x}}(s=',s_pos,')$']);
        hold off

        subplot(1,2,2)
        hold on
        loglog(tplot,(x_act(j,plot_indcs)-x_est(j,plot_indcs)),[colors{j},'o-'],'LineWidth',1.5,'DisplayName',['$\mathbf{e}(s=',s_pos,')$']);
        hold off
    end

    % Clean up the figure
    ax1 = subplot(1,2,1);   ax1.XScale = 'log';     ax1.XTickLabels = {'0.001';'0.01';'0.1';'1'};
    lgd1 = legend('Interpreter','latex');                lgd1.FontSize = 10.5;
    lgd1.Location = 'northwest';
    xlabel('$t$','FontSize',15,'Interpreter','latex');  ylabel('$\mathbf{x}$','FontSize',15,'Interpreter','latex');
    title('PDE state value $\mathbf{x}(t)$ and estimate $\mathbf{\hat{x}}(t)$','Interpreter','latex','FontSize',15);
    ax2 = subplot(1,2,2);   ax2.XScale = 'log';     ax2.XTickLabels = {'0.001';'0.01';'0.1';'1'};
    lgd2 = legend('Interpreter','latex');                lgd2.FontSize = 10.5;
    xlabel('$t$','FontSize',15,'Interpreter','latex');    ylabel('$\mathbf{e}$','FontSize',15,'Interpreter','latex');
    title('Error $\mathbf{e}=\mathbf{\hat{x}}(t)-\mathbf{x}(t)$','Interpreter','latex','FontSize',15);
    fig1.Position = [700 600 1000 450];

    end

    %% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
    function plot_figure_Hinf_optimal_estimator_DEMO(solution,grid,opts,grid_idcs)
    % % % Plot simulated values of actual state and estimated state at several
    % % % grid points, as specified by grid_idcs.

    % Extract actual and estimated solution at each time step.
    x_act = reshape(solution.timedep.pde(:,1,:),opts.N+1,[]);
    x_est = reshape(solution.timedep.pde(:,2,:),opts.N+1,[]);
    tval = solution.timedep.dtime;

    % Set options for the plot
    plot_indcs = floor(logspace(0,log(opts.tf/opts.dt)/log(10),40));
    %plot_indcs = floor(linspace(1,opts.tf/opts.dt,66));
    tplot = tval(plot_indcs);           % Only plot at select times
    colors = {'b','g','m','r','k','c','r','y','o'};     % Colors for the plot

    % Plot evolution of actual and estimated
    fig1 = figure(1);
    for j=grid_idcs
        s_pos = num2str(grid.phys(j));  % Position associated to grid index.
        subplot(1,2,1)
        hold on
        plot(tplot,x_act(j,plot_indcs),[colors{j},'-'],'LineWidth',2,'DisplayName',['$\mathbf{x}(s=',s_pos,')$']);
        plot(tplot,x_est(j,plot_indcs),[colors{j},'o--'],'LineWidth',1.5,'DisplayName',['$\mathbf{\hat{x}}(s=',s_pos,')$']);
        hold off

        subplot(1,2,2)
        hold on
        loglog(tplot,(x_act(j,plot_indcs)-x_est(j,plot_indcs)),[colors{j},'o-'],'LineWidth',1.5,'DisplayName',['$\mathbf{e}(s=',s_pos,')$']);
        hold off
    end
        s_pos = num2str(grid.phys(j));  % Position associated to grid index.
        subplot(1,2,1)
        hold on
        plot(tplot,x_act(j,plot_indcs),[colors{j},'-'],'LineWidth',2,'DisplayName',['$\mathbf{x}(s=',s_pos,')$']);
        plot(tplot,x_est(j,plot_indcs),[colors{j},'o--'],'LineWidth',1.5,'DisplayName',['$\mathbf{\hat{x}}(s=',s_pos,')$']);
        hold off

        subplot(1,2,2)
        hold on
        loglog(tplot,(x_act(j,plot_indcs)-x_est(j,plot_indcs)),[colors{j},'o-'],'LineWidth',1.5,'DisplayName',['$\mathbf{e}(s=',s_pos,')$']);
        hold off
    end
        s_pos = num2str(grid.phys(j));  % Position associated to grid index.
        subplot(1,2,1)
        hold on
        plot(tplot,x_act(j,plot_indcs),[colors{j},'-'],'LineWidth',2,'DisplayName',['$\mathbf{x}(s=',s_pos,')$']);
        plot(tplot,x_est(j,plot_indcs),[colors{j},'o--'],'LineWidth',1.5,'DisplayName',['$\mathbf{\hat{x}}(s=',s_pos,')$']);
        hold off

        subplot(1,2,2)
        hold on
        loglog(tplot,(x_act(j,plot_indcs)-x_est(j,plot_indcs)),[colors{j},'o-'],'LineWidth',1.5,'DisplayName',['$\mathbf{e}(s=',s_pos,')$']);
        hold off
    end

    % Clean up the figure
    ax1 = subplot(1,2,1);   ax1.XScale = 'log';     ax1.XTickLabels = {'0.001';'0.01';'0.1';'1'};
    lgd1 = legend('Interpreter','latex');                lgd1.FontSize = 10.5;
    lgd1.Location = 'northwest';
    xlabel('$t$','FontSize',15,'Interpreter','latex');  ylabel('$\mathbf{x}$','FontSize',15,'Interpreter','latex');
    title('PDE state value $\mathbf{x}(t)$ and estimate $\mathbf{\hat{x}}(t)$','Interpreter','latex','FontSize',15);
    ax2 = subplot(1,2,2);   ax2.XScale = 'log';     ax2.XTickLabels = {'0.001';'0.01';'0.1';'1'};
    lgd2 = legend('Interpreter','latex');                lgd2.FontSize = 10.5;
    xlabel('$t$','FontSize',15,'Interpreter','latex');    ylabel('$\mathbf{e}$','FontSize',15,'Interpreter','latex');
    title('Error $\mathbf{e}=\mathbf{\hat{x}}(t)-\mathbf{x}(t)$','Interpreter','latex','FontSize',15);
    fig1.Position = [700 600 1000 450];

    end