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 (c) 2009-2025 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 = false; 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 = false; % set operation opts = struct(); [oper, opts] = operatormanager(opts, 'default');
global options
opts.norm = 'fro';
RADI parameters
opts.shifts.naive_update_mode = false; opts.radi.compute_sol_fac = true; opts.radi.compute_res = false; opts.radi.get_ZZt = true; 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 = false; 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 mess_err(opts, '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.520281e+03 + 0.000000e+00i normalized residual: 1.095846e-02 relative change in Z: 9.569335e-03 RADI step: 8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.090575e-02 relative change in Z: 1.157454e-02 RADI step: 9 pc: -1.077506e+03 + 0.000000e+00i normalized residual: 1.053820e-03 relative change in Z: 4.404639e-03 RADI step: 10 pc: -4.471525e+00 + 0.000000e+00i normalized residual: 1.040873e-03 relative change in Z: 1.394344e-03 RADI step: 11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 9.980822e-04 relative change in Z: 1.034868e-03 RADI step: 12 pc: -6.526749e+03 + 0.000000e+00i normalized residual: 1.373945e-04 relative change in Z: 1.081194e-03 RADI step: 13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.371014e-04 relative change in Z: 1.115304e-04 RADI step: 14 pc: -1.041276e+03 + 0.000000e+00i normalized residual: 2.386307e-05 relative change in Z: 4.146519e-04 RADI step: 15 pc: -4.471514e+00 + 0.000000e+00i normalized residual: 2.363546e-05 relative change in Z: 1.572055e-04 RADI step: 16 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 2.268074e-05 relative change in Z: 1.171042e-04 RADI step: 18 pc: -6.258103e+03 + -3.947523e+02i normalized residual: 5.405305e-07 relative change in Z: 1.630778e-04 RADI step: 19 pc: -1.117898e+00 + 0.000000e+00i normalized residual: 5.390402e-07 relative change in Z: 1.259984e-05 RADI step: 20 pc: -1.091178e+03 + 0.000000e+00i normalized residual: 7.032812e-08 relative change in Z: 2.683536e-05 RADI step: 21 pc: -1.216261e+03 + 0.000000e+00i normalized residual: 1.149965e-08 relative change in Z: 1.042322e-05 RADI step: 22 pc: -1.006091e+01 + 0.000000e+00i normalized residual: 9.455159e-09 relative change in Z: 2.053063e-05 RADI step: 24 pc: -6.290886e+03 + -1.037961e+03i normalized residual: 4.509840e-10 relative change in Z: 3.297730e-06 RADI step: 25 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 2.371104e-10 relative change in Z: 5.642218e-06 RADI step: 26 pc: -1.222450e+03 + 0.000000e+00i normalized residual: 2.903984e-11 relative change in Z: 5.232219e-07 Size outReg.Z: 905 x 32 Warning: resopts.res takes precedence over resopts; please use resopts as documented ↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a> ↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a> ↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a> ↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/FDM/LQG_FDM_unstable_radi.m',135)">LQG_FDM_unstable_radi (line 135)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/norms/mess_res2_norms.m',113)">mess_res2_norms (line 113)</a> solving the regulator Riccati equation took 0.19 seconds Residual computation -- RADI: 2.903984e-11 | mess_res2_norms: 8.155423e-11 | eigs: 8.155418e-11

Solve filter Riccati equation.
fprintf('Solve filter Riccati equation\n'); fprintf('-----------------------------\n'); eqn.type = 'N'; opts.LDL_T = true; opts.radi.Z0 = ZC0; opts.radi.Y0 = YC0; % Some additional scaling for non-trivial LDL' term. eqn.B = eqn.B * diag(1 ./ sqrt(1:n_unstable)); eqn.R = diag(1:n_unstable); eqn.Q = eye(size(eqn.C, 1)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Set corresponding residual. opts.radi.W0 = WC0 * diag(1 ./ sqrt(1:n_unstable)); opts.radi.T0 = 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.R) * eqn.B'); 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.R * (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 mess_err(opts, '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.014896e+03 + 0.000000e+00i normalized residual: 1.165382e-02 relative change in Z: 1.054083e-02 RADI step: 8 pc: -1.116444e+00 + 0.000000e+00i normalized residual: 1.162231e-02 relative change in Z: 8.254271e-03 RADI step: 10 pc: -6.219560e+02 + -5.812256e+02i normalized residual: 1.580419e-02 relative change in Z: 9.501003e-03 RADI step: 11 pc: -2.794687e+01 + 0.000000e+00i normalized residual: 1.451894e-02 relative change in Z: 3.122320e-03 RADI step: 12 pc: -5.945165e+03 + 0.000000e+00i normalized residual: 1.231246e-02 relative change in Z: 5.577753e-03 RADI step: 13 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 1.228162e-02 relative change in Z: 3.862549e-04 RADI step: 15 pc: -6.293837e+02 + -1.465574e+03i normalized residual: 5.688371e-04 relative change in Z: 5.880379e-03 RADI step: 16 pc: -2.794688e+01 + 0.000000e+00i normalized residual: 5.596881e-04 relative change in Z: 7.449497e-04 RADI step: 18 pc: -9.118909e+02 + -3.038603e+03i normalized residual: 4.287436e-04 relative change in Z: 1.540944e-03 RADI step: 19 pc: -1.117897e+00 + 0.000000e+00i normalized residual: 4.266501e-04 relative change in Z: 1.144259e-04 RADI step: 21 pc: -1.054685e+03 + -1.036088e+02i normalized residual: 6.088963e-05 relative change in Z: 1.274954e-03 RADI step: 22 pc: -1.107356e+03 + 0.000000e+00i normalized residual: 6.080957e-05 relative change in Z: 2.998477e-04 RADI step: 23 pc: -4.472160e+00 + 0.000000e+00i normalized residual: 6.061245e-05 relative change in Z: 7.767297e-05 RADI step: 25 pc: -9.195846e+02 + -3.884575e+03i normalized residual: 2.211912e-05 relative change in Z: 2.110609e-04 RADI step: 26 pc: -1.788599e+01 + 0.000000e+00i normalized residual: 2.193512e-05 relative change in Z: 5.462805e-05 RADI step: 27 pc: -9.428568e+02 + 0.000000e+00i normalized residual: 1.097522e-05 relative change in Z: 1.693259e-04 RADI step: 28 pc: -4.471521e+00 + 0.000000e+00i normalized residual: 1.092157e-05 relative change in Z: 1.433438e-05 RADI step: 30 pc: -8.604715e+02 + -4.586307e+03i normalized residual: 5.185893e-06 relative change in Z: 1.179754e-04 RADI step: 31 pc: -1.788603e+01 + 0.000000e+00i normalized residual: 5.057562e-06 relative change in Z: 1.764861e-05 RADI step: 32 pc: -1.670155e+03 + 0.000000e+00i normalized residual: 1.572012e-06 relative change in Z: 1.212355e-04 RADI step: 33 pc: -9.512802e+00 + 0.000000e+00i normalized residual: 1.573741e-06 relative change in Z: 7.462918e-06 RADI step: 35 pc: -1.020421e+03 + -5.345332e+03i normalized residual: 6.442747e-07 relative change in Z: 4.423560e-05 RADI step: 37 pc: -9.965875e+02 + -3.708754e+03i normalized residual: 4.253481e-07 relative change in Z: 3.217756e-05 RADI step: 38 pc: -2.058356e+03 + 0.000000e+00i normalized residual: 2.126599e-07 relative change in Z: 3.623041e-05 RADI step: 39 pc: -1.006090e+01 + 0.000000e+00i normalized residual: 2.117168e-07 relative change in Z: 1.349507e-06 RADI step: 41 pc: -1.319279e+03 + -6.671383e+03i normalized residual: 1.285788e-07 relative change in Z: 8.884593e-06 RADI step: 43 pc: -1.224416e+03 + -3.118018e+03i normalized residual: 6.424611e-09 relative change in Z: 1.133633e-05 RADI step: 44 pc: -2.626692e+03 + 0.000000e+00i normalized residual: 4.907051e-09 relative change in Z: 3.409770e-06 RADI step: 45 pc: -1.453144e+01 + 0.000000e+00i normalized residual: 4.810975e-09 relative change in Z: 4.329954e-07 RADI step: 47 pc: -1.154669e+03 + -7.261482e+03i normalized residual: 3.292811e-09 relative change in Z: 1.656437e-06 RADI step: 49 pc: -1.068737e+03 + -2.897572e+03i normalized residual: 4.664451e-10 relative change in Z: 2.159986e-06 RADI step: 51 pc: -3.388976e+03 + -4.327969e+02i normalized residual: 1.207816e-10 relative change in Z: 8.067781e-07 RADI step: 53 pc: -1.282176e+03 + -8.496703e+03i normalized residual: 1.921627e-11 relative change in Z: 3.357603e-07 Size outFil.Z: 905 x 73 Warning: resopts.res takes precedence over resopts; please use resopts as documented ↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a> ↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a> ↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a> ↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/FDM/LQG_FDM_unstable_radi.m',193)">LQG_FDM_unstable_radi (line 193)</a> ↳ In <a href="matlab:opentoline('/builds/mess/mmess/norms/mess_res2_norms.m',113)">mess_res2_norms (line 113)</a> solving the filter Riccati equation took 0.32 seconds Residual computations -- RADI: 1.921627e-11 | mess_res2_norms: 1.107844e-10 | eigs: 1.107815e-10
