
% a script to set-up the 1D Ginzburg-Landau equation
% and control it, 
% checking performance using stochastic analysis

%%%% system parameters
u=1; % mean flow
g=1+i*0.1; % diffusion/dispersion
mu0=0.3;  % bifurcation parameter
mu2=-0.01; % parrallelity parameter
N=100; % number of grid points
sss=0.7; % width parameter
xw=-9; % excitation location
xu=0; % actuator location
xy=-4; % sensor location

%%%%  differentiation matrices and Ginzburg-Landau
[x, DM] = herdif(N, 2,0.6); % Hermite differentiation matrices
D1=DM(:,:,1); % first order differentiation
D2=DM(:,:,2); % second order differentiation
mu=mu0+mu2*x.^2/2; % destabilization term
A=-u*D1+diag(mu)+g*D2;% build dynamics

%%%% trapezoidal integration weights
w=([diff(x);0]+[0;diff(x)])/2; % integration weights
Q=diag(w); % inner product matrix

%%%% control operators
fw=exp(-((x-xw)/sss).^2); B1=fw;%excitations
fu=exp(-((x-xu)/sss).^2); B2=fu;% actuation
fy=exp(-((x-xy)/sss).^2);C=(fy.*w).';%sensor

%%%% control parameters
G=0.1;% sensor noise variance
R=50;% control penalty

%%%% compute control and estimation gains
Pc=myric(A,B2*inv(R)*B2',Q); % solve control
K=-inv(R)*B2'*Pc; % extract control gain K
Pe=myric(A',C'*inv(G)*C,B1*B1'); % solve estimation
L=-Pe*C'*inv(G); % extract estimation gain L

%%%% build closed-loop system
AA=[A B2*K;-L*C A+B2*K+L*C]; % closed-loop dynamic matrix
BB=[B1,zeros(N,1);zeros(N,1),-L]; % input for dist and sensor noise

%%%% compute state covariances
H=[1,0;0,G]; % covariance of excitations and sensor noise
P=mylyap(A,A',B1*B1');% state covariance, no control
Pc=mylyap(AA,AA',BB*H*BB');% state covariance, controlled 

%%%% compute rms, and plot results
rms_u=sqrt(real(K*Pc(1:N,1:N)*K')); % actuation rms
rms_nocont=sqrt(real(diag(P))); % no control, state rms
rms_cont=sqrt(real(diag(Pc(1:N,1:N)))); % controlled, state rms

plot(x,fw,'k',x,fy,'m',x,fu*rms_u,'c', ...
  x,mu,'k--',x,rms_nocont,'b',x,rms_cont,'r');
legend('excitations','sensor','actuator', ...
  'mu','rms, no control','rms, control');
