Contents
function LQR_DAE2_SO(istest)
% Computes a Riccati feedback control for a constrained variant of the % triple chain oscillator model from [1] via the Newton-ADI and RADI % methods. The code is using the ideas for second order matrix exploitation % described in [2,3] and second order DAEs from [4]. % % Input: % istest decides whether the function runs as an interactive demo or a % continuous integration test. % (optional; defaults to 0, i.e. interactive demo) % % References % [1] N. Truhar, K. Veselić, An efficient method for estimating the % optimal dampers viscosity for linear vibrating systems using % Lyapunov equation, SIAM J. Matrix Anal. Appl. 31 (1) (2009) 18--39. % https://doi.org/10.1137/070683052 % % [2] P. Benner, J. Saak, Efficient Balancing based MOR for Second Order % Systems Arising in Control of Machine Tools, in: I. Troch, % F. Breitenecker (Eds.), Proceedings of the MathMod 2009, no. 35 in % ARGESIM-Reports, Vienna Univ. of Technology, ARGE Simulation News, % Vienna, Austria, 2009, pp. 1232--1243, https://doi.org/10.11128/arep.35 % % [3] P. Benner, P. Kürschner, J. Saak, Improved second-order balanced % truncation for symmetric systems, IFAC Proceedings Volumes (7th % Vienna International Conference on Mathematical Modelling) 45 (2) % (2012) 758--762. https://doi.org/10.3182/20120215-3-AT-3016.00134 % % [4] M. M. Uddin, Computational methods for model reduction of large-scale % sparse structured descriptor systems, Dissertation, % Otto-von-Guericke-Universität, Magdeburg, Germany (2015). % URL http://nbn-resolving.de/urn:nbn:de:gbv:ma9:1-6535 % % 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) %
if nargin < 1 istest = false; end
set operation
opts = struct();
[oper, opts] = operatormanager(opts, 'dae_2_so');
load problem data
generate problem
nv = 500; np = 10; nin = 2; nout = 3; p = 0.2; alpha = 0.1; beta = 0.1; v = 5; [M, D, K] = triplechain_MSD(nv, alpha, beta, v); nv = size(M, 1); G = zeros(nv, np); while not(rank(full(G)) == np) G = sprand(np, nv, p); end eqn.M_ = M; eqn.E_ = -D; eqn.K_ = -K; eqn.G_ = G; eqn.haveE = true; eqn.alpha = -0.02; eqn.B = [zeros(nv, nin); rand(nv, nin)]; eqn.C = [zeros(nout, nv), rand(nout, nv)]; eqn.type = 'N';
condition numbers of generated input data
mess_fprintf(opts, 'Condition numbers of the generated input data\n'); As = full([zeros(nv, nv), eye(nv, nv), zeros(nv, np); ... K, D, G'; ... G, zeros(np, nv), zeros(np, np)]); mess_fprintf(opts, 'cond(M)=%e\n', condest(eqn.M_)); mess_fprintf(opts, 'cond(D)=%e\n', condest(eqn.E_)); mess_fprintf(opts, 'cond(K)=%e\n', condest(eqn.K_)); mess_fprintf(opts, 'cond(A)=%e\n\n', condest(As));
Condition numbers of the generated input data cond(M)=1.000000e+01 cond(D)=1.220000e+02 cond(K)=3.514842e+06 cond(A)=4.520750e+05
options
opts.norm = 'fro'; % ADI options opts.adi.maxiter = 1000; opts.adi.res_tol = 1e-15; opts.adi.rel_diff_tol = 1e-16; opts.adi.info = 0; opts.shifts.num_desired = 25; opts.shifts.num_Ritz = 50; opts.shifts.num_hRitz = 25; opts.shifts.b0 = ones(2 * nv + np, 1); opts.shifts.method = 'projection'; opts.shifts.num_desired = 6; % Newton options and maximum iteration number opts.nm.maxiter = 20; opts.nm.res_tol = 1e-10; opts.nm.rel_diff_tol = 1e-16; opts.nm.info = 1; opts.nm.accumulateRes = false; opts.nm.linesearch = true; opts.nm.projection.freq = 0; opts.nm.projection.ortho = true; opts.nm.res = struct('maxiter', 10, ... 'tol', 1e-6, ... 'orth', 0);
the actual Newton call
eqn.type = 'T'; t_mess_lrnm = tic; outnm = mess_lrnm(eqn, opts, oper); t_elapsed1 = toc(t_mess_lrnm); mess_fprintf(opts, 'mess_lrnm took %6.2f seconds \n', t_elapsed1); if istest if min(outnm.res) >= opts.nm.res_tol mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result'); end else figure(); semilogy(outnm.res, 'LineWidth', 3); title('0 = C^TC + A^T X E + E^T X A -E^T X BB^T X M'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end [mZ, nZ] = size(outnm.Z); mess_fprintf(opts, 'outnm.Z: %d x%d\n', mZ, nZ);
Using line search (res: 1.541695e+04) lambda: 3.712204e-03 NM step: 1 normalized residual: 9.751872e-01 relative change in K: 1.000000e+00 number of ADI steps: 212 Using line search (res: 2.722341e+01) lambda: 9.427211e-02 NM step: 2 normalized residual: 8.948931e-01 relative change in K: 1.453576e+00 number of ADI steps: 238 Using line search (res: 1.000567e+00) lambda: 5.912759e-01 NM step: 3 normalized residual: 2.666995e-01 relative change in K: 7.611568e-01 number of ADI steps: 258 NM step: 4 normalized residual: 1.644840e-02 relative change in K: 1.452565e-01 number of ADI steps: 243 NM step: 5 normalized residual: 1.477994e-04 relative change in K: 1.187319e-02 number of ADI steps: 240 NM step: 6 normalized residual: 6.317262e-08 relative change in K: 2.433063e-04 number of ADI steps: 256 NM step: 7 normalized residual: 9.542599e-15 relative change in K: 9.455048e-08 number of ADI steps: 276 mess_lrnm took 16.12 seconds outnm.Z: 3002 x333

Lets try the RADI method and compare
opts.norm = 2; % RADI-MESS settings opts.shifts.history = opts.shifts.num_desired * size(eqn.C, 1); opts.shifts.num_desired = 5; % choose either of the three shift methods, here % opts.shifts.method = 'gen-ham-opti'; % opts.shifts.method = 'heur'; opts.shifts.method = 'projection'; opts.shifts.naive_update_mode = false; % .. Suggest false (smart update is faster; convergence is the same). opts.radi.compute_sol_fac = true; opts.radi.maxiter = opts.nm.maxiter * opts.adi.maxiter; opts.radi.get_ZZt = true; opts.radi.res_tol = opts.nm.res_tol; opts.radi.rel_diff_tol = 0; opts.radi.info = 1; t_mess_lrradi = tic; outradi = mess_lrradi(eqn, opts, oper); t_elapsed2 = toc(t_mess_lrradi); mess_fprintf(opts, 'mess_lrradi took %6.2f seconds \n', t_elapsed2); if istest if min(outradi.res) >= opts.radi.res_tol mess_err(opts, 'TEST:accuracy', ... 'unexpectedly inaccurate result'); end else figure(); semilogy(outradi.res, 'LineWidth', 3); title('0 = C^T C + A^T X E + E^T X A -E^T X BB^T X E'); xlabel('number of iterations'); ylabel('normalized residual norm'); end [mZ, nZ] = size(outradi.Z); mess_fprintf(opts, 'outradi.Z: %d x %d\n', mZ, nZ);
RADI step: 1 normalized residual: 9.995782e-01 RADI step: 2 normalized residual: 1.132757e+00 RADI step: 3 normalized residual: 1.193551e+00 RADI step: 4 normalized residual: 1.273884e+00 RADI step: 6 normalized residual: 1.200069e+00 RADI step: 7 normalized residual: 1.110935e+00 RADI step: 8 normalized residual: 9.924135e-01 RADI step: 9 normalized residual: 9.914342e-01 RADI step: 11 normalized residual: 9.859265e-01 RADI step: 13 normalized residual: 9.854293e-01 RADI step: 14 normalized residual: 9.850829e-01 RADI step: 15 normalized residual: 8.989245e-01 RADI step: 17 normalized residual: 8.978241e-01 RADI step: 19 normalized residual: 8.903555e-01 RADI step: 20 normalized residual: 8.871879e-01 RADI step: 21 normalized residual: 3.998905e-01 RADI step: 23 normalized residual: 3.908068e-01 RADI step: 25 normalized residual: 3.875977e-01 RADI step: 26 normalized residual: 3.869476e-01 RADI step: 27 normalized residual: 1.087876e-02 RADI step: 29 normalized residual: 7.627028e-03 RADI step: 31 normalized residual: 7.056055e-03 RADI step: 32 normalized residual: 7.055117e-03 RADI step: 33 normalized residual: 5.971528e-03 RADI step: 35 normalized residual: 4.886775e-03 RADI step: 37 normalized residual: 4.089779e-03 RADI step: 38 normalized residual: 3.680049e-03 RADI step: 39 normalized residual: 3.669507e-03 RADI step: 41 normalized residual: 3.593566e-03 RADI step: 43 normalized residual: 2.951882e-03 RADI step: 45 normalized residual: 4.514321e-04 RADI step: 47 normalized residual: 3.257273e-04 RADI step: 49 normalized residual: 3.204141e-04 RADI step: 51 normalized residual: 3.174169e-04 RADI step: 53 normalized residual: 1.354249e-04 RADI step: 55 normalized residual: 1.007473e-04 RADI step: 56 normalized residual: 9.035894e-05 RADI step: 58 normalized residual: 6.174508e-05 RADI step: 60 normalized residual: 4.345328e-05 RADI step: 61 normalized residual: 4.012420e-05 RADI step: 63 normalized residual: 3.737123e-05 RADI step: 65 normalized residual: 3.287971e-05 RADI step: 66 normalized residual: 3.028791e-05 RADI step: 68 normalized residual: 2.399139e-05 RADI step: 70 normalized residual: 2.270283e-05 RADI step: 71 normalized residual: 2.197402e-05 RADI step: 73 normalized residual: 1.262364e-05 RADI step: 75 normalized residual: 1.206112e-05 RADI step: 76 normalized residual: 1.176018e-05 RADI step: 78 normalized residual: 5.279393e-06 RADI step: 79 normalized residual: 5.276961e-06 RADI step: 81 normalized residual: 4.309727e-06 RADI step: 82 normalized residual: 4.192502e-06 RADI step: 83 normalized residual: 4.186397e-06 RADI step: 85 normalized residual: 2.988854e-06 RADI step: 87 normalized residual: 2.513385e-06 RADI step: 88 normalized residual: 2.499097e-06 RADI step: 90 normalized residual: 1.550546e-06 RADI step: 92 normalized residual: 1.150373e-06 RADI step: 93 normalized residual: 1.100940e-06 RADI step: 94 normalized residual: 1.098193e-06 RADI step: 96 normalized residual: 7.116742e-07 RADI step: 98 normalized residual: 6.202968e-07 RADI step: 99 normalized residual: 6.194809e-07 RADI step: 101 normalized residual: 4.848715e-07 RADI step: 102 normalized residual: 4.832041e-07 RADI step: 104 normalized residual: 4.395502e-07 RADI step: 105 normalized residual: 4.317268e-07 RADI step: 107 normalized residual: 3.981181e-07 RADI step: 109 normalized residual: 3.753182e-07 RADI step: 110 normalized residual: 3.748458e-07 RADI step: 112 normalized residual: 3.024686e-07 RADI step: 114 normalized residual: 2.893857e-07 RADI step: 115 normalized residual: 2.846754e-07 RADI step: 117 normalized residual: 2.662157e-07 RADI step: 119 normalized residual: 6.161359e-08 RADI step: 121 normalized residual: 5.875751e-08 RADI step: 123 normalized residual: 5.380798e-08 RADI step: 125 normalized residual: 2.070259e-08 RADI step: 126 normalized residual: 1.654460e-08 RADI step: 127 normalized residual: 1.653471e-08 RADI step: 129 normalized residual: 9.903800e-09 RADI step: 131 normalized residual: 8.483887e-09 RADI step: 132 normalized residual: 7.728467e-09 RADI step: 133 normalized residual: 7.720468e-09 RADI step: 135 normalized residual: 2.215695e-09 RADI step: 137 normalized residual: 2.414466e-09 RADI step: 138 normalized residual: 2.439949e-09 RADI step: 140 normalized residual: 1.977412e-09 RADI step: 142 normalized residual: 1.936491e-09 RADI step: 144 normalized residual: 1.902959e-09 RADI step: 146 normalized residual: 6.466186e-10 RADI step: 148 normalized residual: 5.933499e-10 RADI step: 150 normalized residual: 5.576693e-10 RADI step: 152 normalized residual: 3.025145e-10 RADI step: 154 normalized residual: 1.545206e-10 RADI step: 155 normalized residual: 1.485310e-10 RADI step: 157 normalized residual: 1.170179e-10 RADI step: 159 normalized residual: 8.168185e-11 mess_lrradi took 1.42 seconds outradi.Z: 3002 x 331

compare
if not(istest) figure(); ls_nm = cumsum([outnm.adi.niter]); ls_radi = 1:outradi.niter; semilogy(ls_nm, outnm.res, 'k--', ... ls_radi, outradi.res, 'b-', ... 'LineWidth', 3); title('0 = C^T C + A^T X E + E^T X A -E^T X BB^T X E'); xlabel('number of solves with A + p * E'); ylabel('normalized residual norm'); legend('LR-NM', 'RADI'); end
