Contents
% DEMO7_observer_based_control.m % See Chapter 11.7 of the manual for a description. % % This document illustrates construction of an observer-based state feedback % controller. We construct the observer and controller separately % and couple them to find a closed loop system which is then simulated for % different IC and disturbances. % % We consider the wave equation PDE defined by: % PDE: \dot{x}(t,s) = (d^2/ds^2) x(t,s) + lam x(t,s) + w(t) + u(t), s in [0,1]; % Outputs: z(t) = [int_{0}^{1} x(t,s) ds + w(t); u(t)]; % y(t0 = x(t,1); % BCs: 0 = x(t,0) = d/ds x(t,1); % % First we convert the above PDE to a PIE of the form: % [T \dot{v}](t,s) = [A v](t,s) + [B1 w](t,s) + [B2 u](t,s); % z(t) = [C1 v](t) + [D11 w](t) + [D12 u](t); % y(t) = [C2 v](t) + [D21 w](t) + [D22 u](t) % % We find the observer gains L by 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)] % % where L = P^{-1}*Z. % Likewise we find a state feedback control u = Kv by solving % the LPI % % min_{gam,P,Z} gam % s.t. P>=0 % [ -gam*I D11 (C1*P+D12*Z)*T' ] % [ D11' -gam*I B1' ] <= 0 % [ T*(C1*P+D12*Z) B1 (A*P+B2*Z)*T'+T*(A*P+B2*Z)'] % % where K = Z*P^{-1}. % % We simulate the open loop and closed loop response of the PDE for various % IC using PIESIM. %
clc; clear; close all; echo on %%%%%%%%%%%%%%%%%% Start Code %%%%%%%%%%%%%%%%%%
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
% % Declare the PDE, and convert it to a PIE. Declare the PDE using command line parser
pvar s t lam = 5; PDE = sys(); x = state('pde'); w = state('in'); y = state('out'); z = state('out', 2); u = state('in'); eqs = [diff(x,t) == diff(x,s,2) + lam*x + s*w + s*u; z == [int(x,s,[0,1]); u]; % change z = int(x,s,[0,1]) y == subs(x,s,1); subs(x,s,0)==0; subs(diff(x,s),s,1)==0]; PDE = addequation(PDE,eqs); PDE = setControl(PDE,u); PDE = setObserve(PDE,y); display_PDE(PDE); % Compute the associated PIE, and extract the operators. PIE = convert(PDE,'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; C2 = PIE.C2; D21 = PIE.D21; D22 = PIE.D22;
%%%%%%%%%%%%%%%%%% Start Code %%%%%%%%%%%%%%%%%% %% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %% % % % Declare the PDE, and convert it to a PIE. % Declare the PDE using command line parser pvar s t lam = 5; PDE = sys(); Initialized sys() object of type "pde" x = state('pde'); w = state('in'); y = state('out'); z = state('out', 2); u = state('in'); eqs = [diff(x,t) == diff(x,s,2) + lam*x + s*w + s*u; z == [int(x,s,[0,1]); u]; % change z = int(x,s,[0,1]) y == subs(x,s,1); subs(x,s,0)==0; subs(diff(x,s),s,1)==0]; PDE = addequation(PDE,eqs); 6 equations were added to sys() object PDE = setControl(PDE,u); 1 inputs were designated as controlled inputs PDE = setObserve(PDE,y); 1 outputs were designated as observed outputs display_PDE(PDE); ∂ₜ x(t,s) = ∂²ₛ x(t,s) + 5 * x(t,s) + C₁₃(s) * w(t) + C₁₄(s) * u(t); y(t) = x(t,1); z(t) = ₀∫¹[C₃₁ * x(t,s)]ds + C₃₂ * u(t); 0 = - x(t,0); 0 = - ∂ₛ x(t,1); Call "PDE.C{i,j}" to see the value of coefficients Cᵢⱼ as in the displayed equations. % Compute the associated PIE, and extract the operators. PIE = convert(PDE,'pie'); --- 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 PIE = PIE.params; T = PIE.T; A = PIE.A; C1 = PIE.C1; B2 = PIE.B2; B1 = PIE.B1; D11 = PIE.D11; D12 = PIE.D12; C2 = PIE.C2; D21 = PIE.D21; D22 = PIE.D22;
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
% % Compute an optimal observer gains and optimal controller gains for the PIE.
% % Use the predefined Hinf estimator executive function. settings = lpisettings('heavy'); [prog_k, Kval, gam_co_val] = PIETOOLS_Hinf_control(PIE, settings); [prog_l, Lval, gam_ob_val] = PIETOOLS_Hinf_estimator(PIE, settings);
%% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %% % % % Compute an optimal observer gains and optimal controller gains for the PIE. % % Use the predefined Hinf estimator executive function. settings = lpisettings('heavy'); [prog_k, Kval, gam_co_val] = PIETOOLS_Hinf_control(PIE, settings); --- 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: 349 120 The coefficient matrix is not full row rank, numerical problems may occur. 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 6 free variables in a quadratic cone eqs m = 120, order n = 34, dim = 351, blocks = 5 nnz(A) = 1927 + 0, nnz(ADA) = 11136, nnz(L) = 5748 it : b*y gap delta rate t/tP* t/tD* feas cg cg prec 0 : 2.43E+01 0.000 1 : 7.72E-01 6.40E+00 0.000 0.2633 0.9000 0.9000 1.21 1 1 1.2E+01 2 : 1.06E+00 1.85E+00 0.000 0.2897 0.9000 0.9000 1.19 1 1 3.5E+00 3 : 1.11E+00 4.98E-01 0.000 0.2686 0.9000 0.9000 0.96 1 1 9.5E-01 4 : 1.14E+00 1.81E-01 0.000 0.3643 0.9000 0.9000 0.76 1 1 4.0E-01 5 : 1.11E+00 6.46E-02 0.000 0.3557 0.9000 0.9000 0.61 1 1 1.8E-01 6 : 1.10E+00 1.82E-02 0.000 0.2812 0.9000 0.9000 0.56 1 1 6.6E-02 7 : 1.12E+00 5.97E-03 0.000 0.3288 0.9000 0.9000 0.35 1 1 3.3E-02 8 : 1.17E+00 2.13E-03 0.000 0.3574 0.9000 0.9000 0.15 1 1 1.9E-02 9 : 1.22E+00 8.58E-04 0.000 0.4020 0.9000 0.9000 0.12 1 1 1.2E-02 10 : 1.26E+00 3.16E-04 0.000 0.3690 0.9000 0.9000 0.20 1 1 6.5E-03 11 : 1.30E+00 1.22E-04 0.000 0.3859 0.9000 0.9000 0.19 1 1 3.9E-03 12 : 1.34E+00 5.01E-05 0.000 0.4103 0.9000 0.9000 0.28 1 1 2.3E-03 13 : 1.37E+00 2.00E-05 0.000 0.3995 0.9000 0.9000 0.20 1 1 1.4E-03 14 : 1.41E+00 6.95E-06 0.000 0.3474 0.9000 0.9000 0.27 1 1 7.1E-04 15 : 1.44E+00 2.42E-06 0.000 0.3474 0.9000 0.9000 0.25 1 1 3.9E-04 16 : 1.46E+00 8.37E-07 0.000 0.3465 0.9000 0.9000 0.28 1 1 2.0E-04 17 : 1.48E+00 2.87E-07 0.000 0.3427 0.9000 0.9000 0.30 1 1 1.0E-04 18 : 1.50E+00 1.06E-07 0.000 0.3698 0.9000 0.9000 0.35 2 2 5.4E-05 19 : 1.51E+00 3.46E-08 0.000 0.3265 0.8797 0.9000 0.37 2 2 2.9E-05 20 : 1.51E+00 1.31E-08 0.000 0.3796 0.9000 0.8798 0.42 2 2 1.5E-05 21 : 1.52E+00 5.73E-09 0.000 0.4356 0.9000 0.9000 0.41 2 2 9.2E-06 22 : 1.52E+00 2.33E-09 0.000 0.4073 0.9000 0.9000 0.48 2 2 4.9E-06 23 : 1.53E+00 1.14E-09 0.000 0.4896 0.9000 0.9000 0.38 3 2 3.3E-06 24 : 1.53E+00 4.25E-10 0.000 0.3718 0.9000 0.9000 0.43 11 7 1.6E-06 25 : 1.53E+00 1.86E-10 0.000 0.4390 0.9000 0.9000 0.40 3 3 9.6E-07 26 : 1.53E+00 6.72E-11 0.000 0.3604 0.9000 0.9000 0.41 3 3 4.8E-07 27 : 1.53E+00 3.21E-11 0.000 0.4782 0.9000 0.9000 0.39 22 24 2.9E-07 Run into numerical problems. iter seconds digits c*x b*y 27 0.5 4.3 1.5338030035e+00 1.5338835672e+00 |Ax-b| = 4.3e-07, [Ay-c]_+ = 3.5E-09, |x|= 1.1e+04, |y|= 1.0e+06 Detailed timing (sec) Pre IPM Post 3.100E-02 2.350E-01 1.500E-02 Max-norms: ||b||=1, ||c|| = 1, Cholesky |add|=4, |skip| = 61, ||L.L|| = 1.76161e+08. Residual norm: 4.3444e-07 iter: 27 feasratio: 0.3918 pinf: 0 dinf: 0 numerr: 1 timing: [0.0310 0.2350 0.0150] wallsec: 0.2810 cpusec: 0.6562 The closed-loop H-infty norm of the given system is upper bounded by: 1.5338 [prog_l, Lval, gam_ob_val] = PIETOOLS_Hinf_estimator(PIE, settings); --- Executing Search for H_infty Optimal Estimator --- - Declaring Positive Storage Operator variable and indefinite Observer 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: 349 114 The coefficient matrix is not full row rank, numerical problems may occur. 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 6 free variables in a quadratic cone eqs m = 114, order n = 34, dim = 351, blocks = 5 nnz(A) = 1825 + 0, nnz(ADA) = 10364, nnz(L) = 5239 it : b*y gap delta rate t/tP* t/tD* feas cg cg prec 0 : 2.43E+01 0.000 1 : 4.99E-01 6.43E+00 0.000 0.2645 0.9000 0.9000 1.31 1 1 1.2E+01 2 : 4.46E-01 2.16E+00 0.000 0.3364 0.9000 0.9000 1.78 1 1 2.9E+00 3 : 4.53E-01 7.52E-01 0.000 0.3476 0.9000 0.9000 0.84 1 1 1.2E+00 4 : 4.52E-01 2.72E-01 0.000 0.3613 0.9000 0.9000 0.45 1 1 6.0E-01 5 : 4.74E-01 8.40E-02 0.000 0.3091 0.9000 0.9000 0.23 1 1 2.9E-01 6 : 4.06E-01 2.30E-02 0.000 0.2735 0.9000 0.9000 0.24 1 1 1.2E-01 7 : 3.24E-01 6.39E-03 0.000 0.2783 0.9000 0.9000 0.28 1 1 5.2E-02 8 : 2.64E-01 1.98E-03 0.000 0.3104 0.9000 0.9000 0.36 1 1 2.3E-02 9 : 2.39E-01 7.79E-04 0.000 0.3929 0.9000 0.9000 0.41 1 1 1.2E-02 10 : 2.27E-01 2.71E-04 0.000 0.3475 0.9000 0.9000 0.40 1 1 5.7E-03 11 : 2.23E-01 1.15E-04 0.000 0.4261 0.9000 0.9000 0.39 1 1 3.2E-03 12 : 2.25E-01 4.34E-05 0.000 0.3759 0.9000 0.9000 0.39 1 1 1.6E-03 13 : 2.28E-01 2.02E-05 0.000 0.4654 0.9000 0.9000 0.39 1 1 9.8E-04 14 : 2.34E-01 7.47E-06 0.000 0.3700 0.9000 0.9000 0.37 1 1 5.0E-04 15 : 2.37E-01 3.49E-06 0.000 0.4666 0.9000 0.9000 0.38 2 2 3.0E-04 16 : 2.42E-01 1.04E-06 0.000 0.2989 0.8692 0.9000 0.36 2 2 1.4E-04 17 : 2.45E-01 4.50E-07 0.000 0.4316 0.9000 0.9000 0.39 2 2 8.1E-05 18 : 2.49E-01 1.61E-07 0.000 0.3575 0.9000 0.9000 0.37 2 2 4.1E-05 19 : 2.51E-01 6.00E-08 0.000 0.3735 0.9000 0.9000 0.42 2 2 2.0E-05 20 : 2.53E-01 2.24E-08 0.000 0.3729 0.9000 0.9000 0.41 2 2 1.0E-05 21 : 2.54E-01 8.25E-09 0.000 0.3687 0.9000 0.9000 0.46 3 3 5.1E-06 22 : 2.55E-01 3.68E-09 0.000 0.4464 0.9000 0.9000 0.44 3 4 3.0E-06 23 : 2.56E-01 1.55E-09 0.000 0.4203 0.9000 0.9000 0.44 4 4 1.7E-06 24 : 2.57E-01 6.98E-10 0.000 0.4510 0.9000 0.9000 0.43 4 6 9.8E-07 Run into numerical problems. iter seconds digits c*x b*y 24 0.3 3.8 2.5645910306e-01 2.5650002154e-01 |Ax-b| = 1.5e-06, [Ay-c]_+ = 1.3E-08, |x|= 1.8e+03, |y|= 3.1e+04 Detailed timing (sec) Pre IPM Post 1.500E-02 1.250E-01 0.000E+00 Max-norms: ||b||=1, ||c|| = 1, Cholesky |add|=2, |skip| = 58, ||L.L|| = 3.59306e+08. Residual norm: 1.4789e-06 iter: 24 feasratio: 0.4341 pinf: 0 dinf: 0 numerr: 1 timing: [0.0150 0.1250 0] wallsec: 0.1400 cpusec: 0.3594 The H-infty gain from disturbance to error in estimated state is upper bounded by: 0.2565
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
% Construct the operators defining the PIE. T_CL = [T, 0*T; 0*T, T]; A_CL = [A, B2*Kval; -Lval*C2, A+Lval*C2]; B_CL = [B1; Lval*D21]; C_CL = [C1, D12*Kval; 0*C1, C1]; D_CL = [D11; 0*D11]; % Declare the PIE. PIE_CL = pie_struct(); PIE_CL.vars = PIE.vars; PIE_CL.dom = PIE.dom; PIE_CL.T = T_CL; PIE_CL.A = A_CL; PIE_CL.B1 = B_CL; PIE_CL.C1 = C_CL; PIE_CL.D11 = D_CL; PIE_CL = initialize(PIE_CL);
%% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %% % Construct the operators defining the PIE. T_CL = [T, 0*T; 0*T, T]; A_CL = [A, B2*Kval; -Lval*C2, A+Lval*C2]; B_CL = [B1; Lval*D21]; C_CL = [C1, D12*Kval; 0*C1, C1]; D_CL = [D11; 0*D11]; % Declare the PIE. PIE_CL = pie_struct(); PIE_CL.vars = PIE.vars; PIE_CL.dom = PIE.dom; PIE_CL.T = T_CL; PIE_CL.A = A_CL; PIE_CL.B1 = B_CL; PIE_CL.C1 = C_CL; PIE_CL.D11 = D_CL; PIE_CL = initialize(PIE_CL);
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %%
% % 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 % 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 = 10; % Simulate up to t = 10; opts.dt = 1e-2; % Use time step of 10^-2 opts.intScheme=1; % Time-step using Backward Differentiation Formula (BDF) ndiff = [0,0,1]; % The PDE state involves 1 second order differentiable state variables % Simulate the solution to the PIE without controller for different IC. uinput.ic.PDE = [-10*sx]; % IC PIE uinput.w = 5*exp(-st); % disturbance [solution_OL,grid] = PIESIM(PIE,opts,uinput,ndiff); % Simulate the solution to the PIE with controller for different IC and disturbance. ndiff = [0,0,2]; uinput.ic.PDE = [-10*sx; 0]; % IC PIE and observed state uinput.w = 5*exp(-st); % disturbance [solution_CL_a,grid] = PIESIM(PIE_CL,opts,uinput,ndiff); uinput.ic.PDE = [sin(sx*pi/2); 0]; % IC PIE uinput.w = 5*sin(pi*st)./(st+eps); % disturbance [solution_CL_b,grid] = PIESIM(PIE_CL,opts,uinput,ndiff);
%% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - %% % % % 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 % 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 = 10; % Simulate up to t = 10; opts.dt = 1e-2; % Use time step of 10^-2 opts.intScheme=1; % Time-step using Backward Differentiation Formula (BDF) ndiff = [0,0,1]; % The PDE state involves 1 second order differentiable state variables % Simulate the solution to the PIE without controller for different IC. uinput.ic.PDE = [-10*sx]; % IC PIE uinput.w = 5*exp(-st); % disturbance [solution_OL,grid] = PIESIM(PIE,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 1.579 or decreasing an order of the scheme (opts.Norder). Value of a regulated output 1 at a final time 10.000000 is 283114491383.5142 Value of a regulated output 2 at a final time 10.000000 is 0.0000 Value of an observed output 1 at a final time 10.000000 is 444715202083.8552 % Simulate the solution to the PIE with controller for different IC and disturbance. ndiff = [0,0,2]; uinput.ic.PDE = [-10*sx; 0]; % IC PIE and observed state uinput.w = 5*exp(-st); % disturbance [solution_CL_a,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 stable for the given problem. Any observed instabilities must be physical. Value of a regulated output 1 at a final time 10.000000 is 0.0000 Value of a regulated output 2 at a final time 10.000000 is -0.0003 Value of a regulated output 3 at a final time 10.000000 is 0.0000 Value of a regulated output 4 at a final time 10.000000 is 0.0000 uinput.ic.PDE = [sin(sx*pi/2); 0]; % IC PIE uinput.w = 5*sin(pi*st)./(st+eps); % disturbance [solution_CL_b,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 stable for the given problem. Any observed instabilities must be physical. Value of a regulated output 1 at a final time 10.000000 is 0.0085 Value of a regulated output 2 at a final time 10.000000 is 0.1240 Value of a regulated output 3 at a final time 10.000000 is -0.0043 Value of a regulated output 4 at a final time 10.000000 is 0.0000
Extract actual solution at each time step.
tval = solution_OL.timedep.dtime; x_OL = reshape(solution_OL.timedep.pde(:,1,:),opts.N+1,[]); x_CL_a = reshape(solution_CL_a.timedep.pde(:,1,:),opts.N+1,[]); hatx_CL_a = reshape(solution_CL_a.timedep.pde(:,2,:),opts.N+1,[]); x_CL_b = reshape(solution_CL_b.timedep.pde(:,1,:),opts.N+1,[]); hatx_CL_b = reshape(solution_CL_b.timedep.pde(:,2,:),opts.N+1,[]);
%% % Extract actual solution at each time step. tval = solution_OL.timedep.dtime; x_OL = reshape(solution_OL.timedep.pde(:,1,:),opts.N+1,[]); x_CL_a = reshape(solution_CL_a.timedep.pde(:,1,:),opts.N+1,[]); hatx_CL_a = reshape(solution_CL_a.timedep.pde(:,2,:),opts.N+1,[]); x_CL_b = reshape(solution_CL_b.timedep.pde(:,1,:),opts.N+1,[]); hatx_CL_b = reshape(solution_CL_b.timedep.pde(:,2,:),opts.N+1,[]);
Set options for the plot
plot_indcs = floor(logspace(0,log(opts.tf/opts.dt)/log(10),5)); %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 grid_idcs = 1:2:9; % Only plot at a few grid points
%% % Set options for the plot plot_indcs = floor(logspace(0,log(opts.tf/opts.dt)/log(10),5)); %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 grid_idcs = 1:2:9; % Only plot at a few grid points
tsteps = 1:50:length(tval); fig1 = figure(1); surf(tval(tsteps),grid.phys,x_OL(:,tsteps,:),'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row'); h=colorbar ; colormap jet box on ylabel(h,'$|\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('$\mathbf{x}(t,s)$','FontSize',15,'Interpreter','latex'); title('Open loop response with $\mathbf{x}(0,s)=\frac{5}{4}(1-s^2)$, $w(t)=5e^{-t}$','Interpreter','latex','FontSize',15); fig2 = figure(2); subplot(2,1,1); surf(tval(tsteps),grid.phys,x_CL_a(:,tsteps,:),'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row'); h=colorbar ; colormap jet box on ylabel(h,'$|\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('$\mathbf{x}(t,s)$','FontSize',15,'Interpreter','latex'); title('Closed loop response with $\mathbf{x}(0,s)=\frac{5}{4}(1-s^2)$, $w(t)=5e^{-t}$','Interpreter','latex','FontSize',15); subplot(2,1,2); surf(tval(tsteps),grid.phys,x_CL_b(:,tsteps,:),'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row'); h=colorbar ; colormap jet box on ylabel(h,'$|\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('$\mathbf{x}(t,s)$','FontSize',15,'Interpreter','latex'); title('Closed loop response with $\mathbf{x}(0,s)=-\frac{4}{\pi^2}\sin(\frac{\pi}{2}s)$, $w(t)=5\frac{\sin(\pi t)}{t}$','Interpreter','latex','FontSize',15);
%% tsteps = 1:50:length(tval); fig1 = figure(1); surf(tval(tsteps),grid.phys,x_OL(:,tsteps,:),'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row'); h=colorbar ; colormap jet box on ylabel(h,'$|\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('$\mathbf{x}(t,s)$','FontSize',15,'Interpreter','latex'); title('Open loop response with $\mathbf{x}(0,s)=\frac{5}{4}(1-s^2)$, $w(t)=5e^{-t}$','Interpreter','latex','FontSize',15); fig2 = figure(2); subplot(2,1,1); surf(tval(tsteps),grid.phys,x_CL_a(:,tsteps,:),'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row'); h=colorbar ; colormap jet box on ylabel(h,'$|\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('$\mathbf{x}(t,s)$','FontSize',15,'Interpreter','latex'); title('Closed loop response with $\mathbf{x}(0,s)=\frac{5}{4}(1-s^2)$, $w(t)=5e^{-t}$','Interpreter','latex','FontSize',15); subplot(2,1,2); surf(tval(tsteps),grid.phys,x_CL_b(:,tsteps,:),'FaceAlpha',0.75,'Linestyle','--','FaceColor','interp','MeshStyle','row'); h=colorbar ; colormap jet box on ylabel(h,'$|\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('$\mathbf{x}(t,s)$','FontSize',15,'Interpreter','latex'); title('Closed loop response with $\mathbf{x}(0,s)=-\frac{4}{\pi^2}\sin(\frac{\pi}{2}s)$, $w(t)=5\frac{\sin(\pi t)}{t}$','Interpreter','latex','FontSize',15);


XX = linspace(0,10,200); w2_tval = subs(5*sin(pi*st)./(st+eps),tval); w1_tval = subs(5*exp(-st),tval); z_quadrature = double(subs(0.5*sx^2-sx,grid.phys)); k_quadrature = double(subs(Kval.Q1,s,grid.phys)); ZZ1 = trapz(z_quadrature,x_CL_a); %int wa ZZ2 = trapz(z_quadrature,x_CL_b);%int wb ZZ3 = trapz(k_quadrature,hatx_CL_a); %u wa ZZ4 = trapz(k_quadrature,hatx_CL_b); %u wb fig3 = figure(3); subplot(1,2,1); hold on; [YY] = spline(tval,ZZ1,XX); plot(XX,YY,'--o','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = 5e^{-t}$']); [YY] = spline(tval,ZZ2,XX); plot(XX,YY,'--x','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = 5\frac{\sin(\pi t)}{t}$']); hold off; subplot(1,2,2); hold on; [YY] = spline(tval,ZZ3,XX); plot(XX,YY,'--o','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = e^{-t}$']); [YY] = spline(tval,ZZ4,XX); plot(XX,YY,'--x','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = 5\frac{\sin(\pi t)}{t}$']); hold off; ax5 = subplot(1,2,1); set(ax5,'XTick',tval(1:900:end)); lgd1 = legend('Interpreter','latex'); lgd1.FontSize = 10.5; lgd1.Location = 'northeast'; xlabel('$t$','FontSize',15,'Interpreter','latex'); ylabel('$z(t)$','FontSize',15,'Interpreter','latex'); title('Closed loop $z(t) = \int_0^1 \mathbf{x}(t,s) ds$','Interpreter','latex','FontSize',15); ax6 = subplot(1,2,2); set(ax6,'XTick',tval(1:900:end)); lgd1 = legend('Interpreter','latex'); lgd1.FontSize = 10.5; lgd1.Location = 'southeast'; xlabel('$t$','FontSize',15,'Interpreter','latex'); ylabel('$u(t)$','FontSize',15,'Interpreter','latex'); title('Control effort $u(t)$','Interpreter','latex','FontSize',15);
%% XX = linspace(0,10,200); w2_tval = subs(5*sin(pi*st)./(st+eps),tval); w1_tval = subs(5*exp(-st),tval); z_quadrature = double(subs(0.5*sx^2-sx,grid.phys)); k_quadrature = double(subs(Kval.Q1,s,grid.phys)); ZZ1 = trapz(z_quadrature,x_CL_a); %int wa ZZ2 = trapz(z_quadrature,x_CL_b);%int wb ZZ3 = trapz(k_quadrature,hatx_CL_a); %u wa ZZ4 = trapz(k_quadrature,hatx_CL_b); %u wb fig3 = figure(3); subplot(1,2,1); hold on; [YY] = spline(tval,ZZ1,XX); plot(XX,YY,'--o','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = 5e^{-t}$']); [YY] = spline(tval,ZZ2,XX); plot(XX,YY,'--x','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = 5\frac{\sin(\pi t)}{t}$']); hold off; subplot(1,2,2); hold on; [YY] = spline(tval,ZZ3,XX); plot(XX,YY,'--o','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = e^{-t}$']); [YY] = spline(tval,ZZ4,XX); plot(XX,YY,'--x','MarkerIndices',1:90:length(XX),'LineWidth',2,'DisplayName',['$w(t) = 5\frac{\sin(\pi t)}{t}$']); hold off; ax5 = subplot(1,2,1); set(ax5,'XTick',tval(1:900:end)); lgd1 = legend('Interpreter','latex'); lgd1.FontSize = 10.5; lgd1.Location = 'northeast'; xlabel('$t$','FontSize',15,'Interpreter','latex'); ylabel('$z(t)$','FontSize',15,'Interpreter','latex'); title('Closed loop $z(t) = \int_0^1 \mathbf{x}(t,s) ds$','Interpreter','latex','FontSize',15); ax6 = subplot(1,2,2); set(ax6,'XTick',tval(1:900:end)); lgd1 = legend('Interpreter','latex'); lgd1.FontSize = 10.5; lgd1.Location = 'southeast'; xlabel('$t$','FontSize',15,'Interpreter','latex'); ylabel('$u(t)$','FontSize',15,'Interpreter','latex'); title('Control effort $u(t)$','Interpreter','latex','FontSize',15);

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