Contents
function LQG_FDM_unstable_nwt(n0, n_unstable, istest)
% Computes stabilizing and detecting solutions of the regulator and filter % Riccati equations via the Newton-ADI method, 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_unstable, n_stable), eye(n_unstable)]; KB0 = Bt' * Y * Vl; if exist('lyap', 'file') % Solve Lyapunov equation. Au = lyap(Y, -Bt * Bt'); else Au = lyap2solve(Y, -Bt * Bt'); end KC0 = C(:, n_stable + 1:n) * Y * Vl; % 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';
ADI tolerances and maximum iteration number
opts.adi.maxiter = 200; opts.adi.res_tol = 1e-10; opts.adi.rel_diff_tol = 1e-16; opts.adi.info = 0; opts.adi.accumulateDeltaK = true; opts.adi.accumulateK = false; opts.adi.compute_sol_fac = true;
shift parameters via projection
opts.shifts.num_desired = 5;
opts.shifts.method = 'projection';
Newton tolerances and maximum iteration number
opts.nm.maxiter = 15; opts.nm.res_tol = 1e-8; opts.nm.rel_diff_tol = 1e-16; opts.nm.info = 1; opts.nm.linesearch = true; opts.nm.accumulateRes = true; opts.nm.res = struct('maxiter', 10, 'tol', 1e-6, 'orth', 0);
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.nm.K0 = KB0; time_nwt_reguator = tic; [outReg, eqn, opts, oper] = mess_lrnm(eqn, opts, oper); t_elapsed1 = toc(time_nwt_reguator); % 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.nm, []) / 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 -- Newton: %e | ' ... 'mess_res2_norms: %e | eigs: %e \n'], ... outReg.res(end), res1, res2); % Print convergence behavior. if istest if min(outReg.res) >= opts.nm.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 BB^T X'); xlabel('number of iteration steps'); ylabel('normalized residual norm'); pause(1); end fprintf('\n');
Solve regulator Riccati equation -------------------------------- NM step: 1 normalized residual: 9.444737e-01 relative change in K: 4.999150e-01 number of ADI steps: 20 NM step: 2 normalized residual: 1.049167e-01 relative change in K: 1.999310e-01 number of ADI steps: 21 NM step: 3 normalized residual: 2.139893e-03 relative change in K: 2.939241e-02 number of ADI steps: 20 NM step: 4 normalized residual: 9.675243e-07 relative change in K: 6.253776e-04 number of ADI steps: 25 NM step: 5 normalized residual: 1.057959e-09 relative change in K: 2.830132e-07 number of ADI steps: 19 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_nwt.m',130)">LQG_FDM_unstable_nwt (line 130)</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.63 seconds Residual computation -- Newton: 1.057959e-09 | mess_res2_norms: 8.313213e-10 | eigs: 8.313213e-10

Solve filter Riccati equation.
fprintf('Solve filter Riccati equation\n'); fprintf('-----------------------------\n'); eqn.type = 'N'; opts.LDL_T = true; opts.nm.K0 = KC0; % Some additional scaling for non-trivial LDL' term. eqn.T = diag(1:n_unstable); eqn.B = eqn.B * diag(1 ./ sqrt(1:n_unstable)); time_nwt_filter = tic; [outFil, eqn, opts, oper] = mess_lrnm(eqn, opts, oper); t_elapsed2 = toc(time_nwt_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.T) * eqn.B'); res3 = mess_res2_norms(outFil.Z, 'riccati', ... eqn, opts, oper, opts.nm, 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.T * (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 -- Newton: %e | ' ... 'mess_res2_norms: %e | eigs: %e \n'], ... outFil.res(end), res3, res4); % Print convergence behavior. if istest if min(outFil.res) >= opts.nm.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 ----------------------------- NM step: 1 normalized residual: 2.741011e-01 relative change in K: 2.356433e-01 number of ADI steps: 61 NM step: 2 normalized residual: 1.691022e-02 relative change in K: 5.632558e-02 number of ADI steps: 49 NM step: 3 normalized residual: 2.806716e-04 relative change in K: 5.727899e-03 number of ADI steps: 61 NM step: 4 normalized residual: 1.256115e-07 relative change in K: 1.122036e-04 number of ADI steps: 70 NM step: 5 normalized residual: 6.005360e-09 relative change in K: 5.056964e-08 number of ADI steps: 43 Size outFil.Z: 905 x 516 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_nwt.m',182)">LQG_FDM_unstable_nwt (line 182)</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 1.60 seconds Residual computations -- Newton: 6.005360e-09 | mess_res2_norms: 1.856792e-09 | eigs: 8.000000e-01
