Contents
function LQG_FDM_unstable_radi(n0, n_unstable, istest)
% Computes stabilizing and detecting solutions of the regulator and filter % Riccati equations via the RADI, respectively. The solution of the % regulator equation is approximated via a ZZ' factorization, and for the % filter equation an LDL' approximation is used. % Note: The LDL' case is only used for demonstration, but with no % underlying necessity here. % % Inputs: % % n0 n0^2 gives the dimension of the original model, i.e. n0 is % the number of degrees of freedom, i.e. grid points, per % spatial direction % (optional; defaults to 30) % % n_unstable number of unstable eigenvalues by construction % (optional; defaults to 5) % % istest flag to determine whether this demo runs as a CI test or % interactive demo % (optional, defaults to 0, i.e. interactive demo) % % % This file is part of the M-M.E.S.S. project % (http://www.mpi-magdeburg.mpg.de/projects/mess). % Copyright © 2009-2022 Jens Saak, Martin Koehler, Peter Benner and others. % All rights reserved. % License: BSD 2-Clause License (see COPYING) % narginchk(0,3); if nargin<1, n0 = 30; end if nargin<2, n_unstable = 5; end if nargin<3, istest = 0; end
construct unstable matrix
Generate the stable part.
A0 = fdm_2d_matrix(n0, '10*x', '1000*y', '0'); n_stable = size(A0, 1); n = n_stable + n_unstable; B = [eye(n_stable, n_unstable); diag(1:n_unstable)]; C = [transpose(B); zeros(2, n)]; C(n_unstable + 1, n_unstable + 1) = 1; C(n_unstable + 2, n_unstable + 2) = 1; % The instability is added by fixing the solution Y of the partial % stabilization problem via Bernoulli and solving the resulting Lyapunov % equation for the unstable matrix Au: % % Au'*Y + Y*Au - Y*Bt*Bt'*Y = 0. % % The resulting stabilizing feedback is defined such that the unstable % eigenvalues of the system are mirrored into the left open half-plane. Y = eye(n_unstable); Bt = B(n_stable+1:n, : ); % Part of B acting on unstable part. Vl = [zeros(n_stable, n_unstable); eye(n_unstable)]; if exist('lyap','file') % Solve Lyapunov equation. Au = lyap(Y, -Bt * Bt'); else Au = lyap2solve(Y, (-Bt * Bt')'); end ZB0 = Vl; YB0 = Y \ eye(n_unstable); WB0 = C'; ZC0 = ZB0; YC0 = YB0; WC0 = B; % The full A is constructed via additive decomposition (block-diagonal). eqn.A_ = blkdiag(A0,Au); eqn.B = B; eqn.C = C; eqn.haveE = 0; % set operation oper = operatormanager('default');
global options
opts.norm = 'fro';
RADI parameters
opts.shifts.naive_update_mode = 0; opts.radi.compute_sol_fac = 1; opts.radi.compute_res = 0; opts.radi.get_ZZt = 1; opts.radi.maxiter = 200; opts.radi.res_tol = 1.0e-10; opts.radi.rel_diff_tol = 1.0e-16; opts.radi.info = 1; % Only for mess_res2_norms (not needed in RADI): opts.radi.res = struct('maxiter', 10, 'tol', 1.0e-06, 'orth', 0);
shift parameters via projection
opts.shifts.num_desired = 5;
opts.shifts.method = 'projection';
for shift banned eigenvalues
opts.shifts.banned = -eig(Au); opts.shifts.banned_tol = 1e-6;
Solve regulator Riccati equation.
fprintf('Solve regulator Riccati equation\n'); fprintf('--------------------------------\n'); eqn.type = 'T'; opts.LDL_T = 0; opts.radi.Z0 = ZB0; opts.radi.Y0 = YB0; opts.radi.W0 = WB0; time_radi_regulator = tic; [outReg, eqn, opts, oper] = mess_lrradi(eqn, opts, oper); t_elapsed1 = toc(time_radi_regulator); % Size of solution factor. fprintf('\nSize outReg.Z: %d x %d\n', ... size(outReg.Z, 1), size(outReg.Z, 2)); % Check residuals. res0 = norm(eqn.C * eqn.C'); res1 = mess_res2_norms(outReg.Z, 'riccati', ... eqn, opts, oper, opts.radi, []) / res0; res2 = abs(eigs(@(x) eqn.A_' * (outReg.Z * ((outReg.Z' * x))) ... + (outReg.Z * ((outReg.Z' * (eqn.A_ * x)))) ... + eqn.C' * (eqn.C * x) - outReg.K' * (outReg.K * x), ... n, 1, 'LM')) / res0; fprintf('solving the regulator Riccati equation took %6.2f seconds \n', ... t_elapsed1); fprintf(['Residual computation -- RADI: %e | ' ... 'mess_res2_norms: %e | eigs: %e \n'], ... outReg.res(end), res1, res2); % Print convergence behavior. if istest if min(outReg.res) >= opts.radi.res_tol error('MESS:TEST:accuracy', 'unexpectedly inaccurate result'); end else figure(1); semilogy(outReg.res,'LineWidth', 3); title('0 = C^T C + A^T X + X A - X B B^T X'); xlabel('number of iteration steps'); ylabel('normalized residual norm'); pause(1); end fprintf('\n');
Solve regulator Riccati equation -------------------------------- RADI step: 1 pc: -1.250003e+01 + 0.000000e+00i normalized residual: 2.504853e-01 relative change in Z: 6.334954e-01 RADI step: 2 pc: -8.000025e+00 + 0.000000e+00i normalized residual: 1.248903e-01 relative change in Z: 2.307932e-01 RADI step: 3 pc: -4.500025e+00 + 0.000000e+00i normalized residual: 9.837321e-02 relative change in Z: 1.554579e-01 RADI step: 4 pc: -2.000025e+00 + 0.000000e+00i normalized residual: 9.051185e-02 relative change in Z: 1.122526e-01 RADI step: 5 pc: -5.000250e-01 + 0.000000e+00i normalized residual: 8.921764e-02 relative change in Z: 3.235217e-02 RADI step: 6 pc: -2.759799e+01 + 0.000000e+00i normalized residual: 7.456561e-02 relative change in Z: 4.044819e-02 RADI step: 7 pc: -6.494748e+03 + 0.000000e+00i normalized residual: 1.089110e-02 relative change in Z: 9.575758e-03 RADI step: 8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.083860e-02 relative change in Z: 1.157451e-02 RADI step: 9 pc: -1.077506e+03 + 0.000000e+00i normalized residual: 1.040756e-03 relative change in Z: 4.392728e-03 RADI step: 10 pc: -4.471525e+00 + 0.000000e+00i normalized residual: 1.027884e-03 relative change in Z: 1.394280e-03 RADI step: 11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 9.855012e-04 relative change in Z: 1.034566e-03 RADI step: 12 pc: -6.526616e+03 + 0.000000e+00i normalized residual: 1.362372e-04 relative change in Z: 1.074209e-03 RADI step: 13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.359461e-04 relative change in Z: 1.115155e-04 RADI step: 14 pc: -1.041276e+03 + 0.000000e+00i normalized residual: 2.361170e-05 relative change in Z: 4.130731e-04 RADI step: 15 pc: -4.471514e+00 + 0.000000e+00i normalized residual: 2.338554e-05 relative change in Z: 1.571960e-04 RADI step: 16 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 2.243961e-05 relative change in Z: 1.170509e-04 RADI step: 18 pc: -6.227725e+03 + -3.597591e+02i normalized residual: 5.256882e-07 relative change in Z: 1.623687e-04 RADI step: 19 pc: -1.117898e+00 + 0.000000e+00i normalized residual: 5.242264e-07 relative change in Z: 1.259725e-05 RADI step: 20 pc: -1.091629e+03 + 0.000000e+00i normalized residual: 6.797841e-08 relative change in Z: 2.649700e-05 RADI step: 21 pc: -1.216101e+03 + 0.000000e+00i normalized residual: 1.115192e-08 relative change in Z: 1.029177e-05 RADI step: 22 pc: -1.006091e+01 + 0.000000e+00i normalized residual: 9.091590e-09 relative change in Z: 2.052900e-05 RADI step: 24 pc: -6.274322e+03 + -1.064010e+03i normalized residual: 4.449329e-10 relative change in Z: 3.234767e-06 RADI step: 25 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 2.292637e-10 relative change in Z: 5.641699e-06 RADI step: 26 pc: -1.222273e+03 + 0.000000e+00i normalized residual: 2.796066e-11 relative change in Z: 5.144872e-07 Size outReg.Z: 905 x 32 solving the regulator Riccati equation took 0.18 seconds Residual computation -- RADI: 2.796066e-11 | mess_res2_norms: 8.155300e-11 | eigs: 8.155296e-11

Solve filter Riccati equation.
fprintf('Solve filter Riccati equation\n'); fprintf('-----------------------------\n'); eqn.type = 'N'; opts.LDL_T = 1; opts.radi.Z0 = ZC0; opts.radi.Y0 = YC0; % Some additional scaling for non-trivial LDL' term. eqn.S = diag(1:n_unstable); eqn.B = eqn.B * diag(1 ./ sqrt(1:n_unstable)); % Set corresponding residual. opts.radi.W0 = WC0 * diag(1 ./ sqrt(1:n_unstable)); opts.radi.S0 = diag(1:n_unstable); time_radi_filter = tic; [outFil, eqn, opts, oper] = mess_lrradi(eqn, opts, oper); t_elapsed2 = toc(time_radi_filter); % Size of solution factor. fprintf('\nSize outFil.Z: %d x %d\n', ... size(outFil.Z, 1), size(outFil.Z, 2)); % Check residuals. res0 = norm((eqn.B * eqn.S) * eqn.B'); eqn.S_diag = diag(eqn.S); res3 = mess_res2_norms(outFil.Z, 'riccati', ... eqn, opts, oper, opts.radi, outFil.D) / res0; ZD = outFil.Z * outFil.D; res4 = abs(eigs(@(x) eqn.A_ * (ZD * ((outFil.Z' * x))) ... + (ZD * ((outFil.Z' * (eqn.A_' * x)))) ... + eqn.B * (eqn.S * (eqn.B' * x)) - (outFil.K)' * ((outFil.K) * x), ... n, 1, 'LM')) / res0; fprintf('solving the filter Riccati equation took %6.2f seconds \n', ... t_elapsed2); fprintf(['Residual computations -- RADI: %e | ' ... 'mess_res2_norms: %e | eigs: %e \n'], ... outFil.res(end), res3, res4); % Print convergence behavior. if istest if min(outFil.res) >= opts.radi.res_tol error('MESS:TEST:accuracy', 'unexpectedly inaccurate result'); end else figure(2); semilogy(outFil.res,'LineWidth', 3); title('0 = B S B^T + A Y + Y A^T - Y C^T C Y'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end
Solve filter Riccati equation ----------------------------- RADI step: 1 pc: -1.250003e+01 + 0.000000e+00i normalized residual: 2.471852e-01 relative change in Z: 8.425445e-01 RADI step: 2 pc: -8.000025e+00 + 0.000000e+00i normalized residual: 1.179628e-01 relative change in Z: 2.442724e-01 RADI step: 3 pc: -4.500025e+00 + 0.000000e+00i normalized residual: 8.950024e-02 relative change in Z: 1.303531e-01 RADI step: 4 pc: -2.000025e+00 + 0.000000e+00i normalized residual: 8.083576e-02 relative change in Z: 8.548480e-02 RADI step: 5 pc: -5.000250e-01 + 0.000000e+00i normalized residual: 7.940022e-02 relative change in Z: 2.739268e-02 RADI step: 6 pc: -2.759799e+01 + 0.000000e+00i normalized residual: 6.388342e-02 relative change in Z: 6.269812e-02 RADI step: 7 pc: -5.014781e+03 + 0.000000e+00i normalized residual: 1.165367e-02 relative change in Z: 1.054087e-02 RADI step: 8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.162216e-02 relative change in Z: 8.254271e-03 RADI step: 10 pc: -6.220259e+02 + -5.816253e+02i normalized residual: 1.580948e-02 relative change in Z: 9.500773e-03 RADI step: 11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 1.452220e-02 relative change in Z: 3.124170e-03 RADI step: 12 pc: -5.945117e+03 + 0.000000e+00i normalized residual: 1.231813e-02 relative change in Z: 5.578680e-03 RADI step: 13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.228722e-02 relative change in Z: 3.865830e-04 RADI step: 15 pc: -6.293714e+02 + -1.465612e+03i normalized residual: 5.689496e-04 relative change in Z: 5.881489e-03 RADI step: 16 pc: -2.794688e+01 + 0.000000e+00i normalized residual: 5.596639e-04 relative change in Z: 7.459946e-04 RADI step: 18 pc: -9.118950e+02 + -3.038626e+03i normalized residual: 4.300191e-04 relative change in Z: 1.540990e-03 RADI step: 19 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 4.279174e-04 relative change in Z: 1.146041e-04 RADI step: 21 pc: -1.054695e+03 + -1.036002e+02i normalized residual: 6.087122e-05 relative change in Z: 1.277466e-03 RADI step: 22 pc: -1.107375e+03 + 0.000000e+00i normalized residual: 6.078395e-05 relative change in Z: 2.997363e-04 RADI step: 23 pc: -4.472160e+00 + 0.000000e+00i normalized residual: 6.058693e-05 relative change in Z: 7.767462e-05 RADI step: 25 pc: -9.195671e+02 + -3.884581e+03i normalized residual: 2.210906e-05 relative change in Z: 2.110296e-04 RADI step: 26 pc: -1.788599e+01 + 0.000000e+00i normalized residual: 2.192515e-05 relative change in Z: 5.462723e-05 RADI step: 27 pc: -9.430082e+02 + 0.000000e+00i normalized residual: 1.096999e-05 relative change in Z: 1.693018e-04 RADI step: 28 pc: -4.471521e+00 + 0.000000e+00i normalized residual: 1.091638e-05 relative change in Z: 1.433180e-05 RADI step: 30 pc: -8.604730e+02 + -4.586325e+03i normalized residual: 5.182548e-06 relative change in Z: 1.179524e-04 RADI step: 31 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 5.054292e-06 relative change in Z: 1.764528e-05 RADI step: 32 pc: -1.669957e+03 + 0.000000e+00i normalized residual: 1.571658e-06 relative change in Z: 1.212049e-04 RADI step: 33 pc: -9.514423e+00 + 0.000000e+00i normalized residual: 1.573388e-06 relative change in Z: 7.461725e-06 RADI step: 35 pc: -1.020407e+03 + -5.345356e+03i normalized residual: 6.440315e-07 relative change in Z: 4.423136e-05 RADI step: 37 pc: -9.965652e+02 + -3.708808e+03i normalized residual: 4.251783e-07 relative change in Z: 3.217101e-05 RADI step: 38 pc: -2.058298e+03 + 0.000000e+00i normalized residual: 2.126008e-07 relative change in Z: 3.622314e-05 RADI step: 39 pc: -1.006090e+01 + 0.000000e+00i normalized residual: 2.116580e-07 relative change in Z: 1.349319e-06 RADI step: 41 pc: -1.319288e+03 + -6.671412e+03i normalized residual: 1.285451e-07 relative change in Z: 8.883680e-06 RADI step: 43 pc: -1.224419e+03 + -3.118088e+03i normalized residual: 6.423242e-09 relative change in Z: 1.133498e-05 RADI step: 44 pc: -2.626777e+03 + 0.000000e+00i normalized residual: 4.906290e-09 relative change in Z: 3.409488e-06 RADI step: 45 pc: -1.455835e+01 + 0.000000e+00i normalized residual: 4.810056e-09 relative change in Z: 4.333556e-07 RADI step: 47 pc: -1.154675e+03 + -7.261498e+03i normalized residual: 3.292168e-09 relative change in Z: 1.656312e-06 RADI step: 49 pc: -1.068699e+03 + 2.897627e+03i normalized residual: 4.663999e-10 relative change in Z: 2.159784e-06 RADI step: 51 pc: -3.389061e+03 + -4.327642e+02i normalized residual: 1.207785e-10 relative change in Z: 8.067382e-07 RADI step: 53 pc: -1.282176e+03 + -8.496731e+03i normalized residual: 1.921659e-11 relative change in Z: 3.357523e-07 Size outFil.Z: 905 x 73 solving the filter Riccati equation took 0.34 seconds Residual computations -- RADI: 1.921659e-11 | mess_res2_norms: 1.107772e-10 | eigs: 1.107844e-10
