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 © 2009-2022 Jens Saak, Martin Koehler, Peter Benner and others. % All rights reserved. % License: BSD 2-Clause License (see COPYING) %
if nargin<1, istest=0; end
set operation
oper = operatormanager('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 rank(full(G))~=np, G = sprand(np,nv,p);end eqn.M_=M; eqn.E_=-D; eqn.K_=-K; eqn.G_=G; eqn.haveE=1; 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
fprintf('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)]); fprintf('cond(M)=%e\n',condest(eqn.M_)); fprintf('cond(D)=%e\n',condest(eqn.E_)); fprintf('cond(K)=%e\n',condest(eqn.K_)); fprintf('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)=7.564144e+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 = 0; opts.nm.linesearch = 1; opts.nm.projection.freq=0; opts.nm.projection.ortho=1; 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); fprintf(1,'mess_lrnm took %6.2f seconds ',t_elapsed1); if istest if min(outnm.res)>=opts.nm.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); 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 disp('outnm.Z:'); disp(size(outnm.Z));
Using line search (res: 1.571023e+04) lambda: 3.526910e-03 NM step: 1 normalized residual: 9.788537e-01 relative change in K: 1.000000e+00 number of ADI steps: 217 Using line search (res: 3.607387e+01) lambda: 7.817125e-02 NM step: 2 normalized residual: 9.129159e-01 relative change in K: 1.565890e+00 number of ADI steps: 242 Using line search (res: 1.295289e+00) lambda: 5.463382e-01 NM step: 3 normalized residual: 2.979520e-01 relative change in K: 7.881720e-01 number of ADI steps: 269 NM step: 4 normalized residual: 2.340302e-02 relative change in K: 1.726727e-01 number of ADI steps: 279 NM step: 5 normalized residual: 1.296103e-04 relative change in K: 1.129738e-02 number of ADI steps: 249 NM step: 6 normalized residual: 5.954607e-07 relative change in K: 7.367418e-04 number of ADI steps: 252 NM step: 7 normalized residual: 2.161913e-11 relative change in K: 4.443217e-06 number of ADI steps: 246 mess_lrnm took 18.38 seconds outnm.Z: 3002 334

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 = 1; opts.radi.maxiter = opts.nm.maxiter * opts.adi.maxiter; opts.radi.get_ZZt = 1; 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); fprintf(1,'mess_lrradi took %6.2f seconds ' ,t_elapsed2); if istest if min(outradi.res)>=opts.radi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(2); 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 disp('outradi.Z:'); disp(size(outradi.Z));
RADI step: 1 normalized residual: 4.718555e-01 RADI step: 2 normalized residual: 5.177764e-01 RADI step: 3 normalized residual: 5.903238e-01 RADI step: 5 normalized residual: 4.717140e-01 RADI step: 7 normalized residual: 4.704211e-01 RADI step: 8 normalized residual: 4.681679e-01 RADI step: 9 normalized residual: 4.679178e-01 RADI step: 11 normalized residual: 4.659111e-01 RADI step: 13 normalized residual: 4.656046e-01 RADI step: 14 normalized residual: 4.654265e-01 RADI step: 16 normalized residual: 4.635360e-01 RADI step: 18 normalized residual: 4.632129e-01 RADI step: 19 normalized residual: 4.628559e-01 RADI step: 21 normalized residual: 4.628494e-01 RADI step: 23 normalized residual: 4.614028e-01 RADI step: 24 normalized residual: 4.614021e-01 RADI step: 25 normalized residual: 3.939312e-01 RADI step: 27 normalized residual: 3.931306e-01 RADI step: 29 normalized residual: 3.876095e-01 RADI step: 30 normalized residual: 3.876023e-01 RADI step: 31 normalized residual: 1.111521e-01 RADI step: 33 normalized residual: 1.108436e-01 RADI step: 35 normalized residual: 1.091657e-01 RADI step: 36 normalized residual: 1.084699e-01 RADI step: 37 normalized residual: 1.108157e-03 RADI step: 38 normalized residual: 1.105713e-03 RADI step: 40 normalized residual: 1.026253e-03 RADI step: 41 normalized residual: 7.535701e-04 RADI step: 42 normalized residual: 7.523806e-04 RADI step: 43 normalized residual: 7.552420e-04 RADI step: 45 normalized residual: 5.593115e-04 RADI step: 46 normalized residual: 5.591406e-04 RADI step: 47 normalized residual: 5.623375e-04 RADI step: 49 normalized residual: 4.878967e-04 RADI step: 51 normalized residual: 4.446146e-04 RADI step: 52 normalized residual: 3.876765e-04 RADI step: 53 normalized residual: 3.869525e-04 RADI step: 55 normalized residual: 2.639868e-04 RADI step: 57 normalized residual: 2.240680e-04 RADI step: 58 normalized residual: 1.972749e-04 RADI step: 60 normalized residual: 1.048241e-04 RADI step: 61 normalized residual: 1.041285e-04 RADI step: 63 normalized residual: 9.151939e-05 RADI step: 64 normalized residual: 7.626034e-05 RADI step: 66 normalized residual: 5.088659e-05 RADI step: 67 normalized residual: 5.036082e-05 RADI step: 69 normalized residual: 4.138188e-05 RADI step: 70 normalized residual: 3.808337e-05 RADI step: 72 normalized residual: 3.254392e-05 RADI step: 74 normalized residual: 2.777236e-05 RADI step: 75 normalized residual: 2.840423e-05 RADI step: 77 normalized residual: 2.842105e-05 RADI step: 79 normalized residual: 2.339757e-05 RADI step: 80 normalized residual: 2.399427e-05 RADI step: 82 normalized residual: 2.409424e-05 RADI step: 84 normalized residual: 2.509118e-05 RADI step: 86 normalized residual: 2.609361e-05 RADI step: 88 normalized residual: 2.548160e-05 RADI step: 90 normalized residual: 2.736158e-05 RADI step: 91 normalized residual: 3.156695e-05 RADI step: 93 normalized residual: 3.279585e-05 RADI step: 95 normalized residual: 3.294197e-05 RADI step: 96 normalized residual: 3.647899e-05 RADI step: 98 normalized residual: 3.603516e-05 RADI step: 100 normalized residual: 3.605483e-05 RADI step: 101 normalized residual: 3.605722e-05 RADI step: 103 normalized residual: 3.565431e-05 RADI step: 104 normalized residual: 3.564548e-05 RADI step: 106 normalized residual: 3.551901e-05 RADI step: 107 normalized residual: 3.453377e-05 RADI step: 108 normalized residual: 3.453358e-05 RADI step: 110 normalized residual: 3.445112e-05 RADI step: 112 normalized residual: 3.312511e-05 RADI step: 113 normalized residual: 3.244876e-05 RADI step: 114 normalized residual: 3.244829e-05 RADI step: 116 normalized residual: 2.429289e-05 RADI step: 118 normalized residual: 2.286353e-05 RADI step: 119 normalized residual: 2.181939e-05 RADI step: 120 normalized residual: 2.181794e-05 RADI step: 122 normalized residual: 1.407987e-05 RADI step: 124 normalized residual: 1.361069e-05 RADI step: 125 normalized residual: 1.333658e-05 RADI step: 127 normalized residual: 4.173386e-06 RADI step: 128 normalized residual: 4.171614e-06 RADI step: 130 normalized residual: 3.863868e-06 RADI step: 131 normalized residual: 3.804309e-06 RADI step: 133 normalized residual: 8.701859e-08 RADI step: 135 normalized residual: 7.061379e-08 RADI step: 136 normalized residual: 6.769437e-08 RADI step: 138 normalized residual: 2.322745e-08 RADI step: 140 normalized residual: 1.984084e-08 RADI step: 142 normalized residual: 1.862937e-08 RADI step: 144 normalized residual: 1.851489e-08 RADI step: 146 normalized residual: 1.785414e-08 RADI step: 147 normalized residual: 1.579542e-08 RADI step: 149 normalized residual: 1.539676e-08 RADI step: 151 normalized residual: 1.382096e-08 RADI step: 153 normalized residual: 1.157142e-08 RADI step: 155 normalized residual: 9.791703e-09 RADI step: 157 normalized residual: 9.744947e-09 RADI step: 158 normalized residual: 9.626649e-09 RADI step: 160 normalized residual: 8.646148e-09 RADI step: 161 normalized residual: 8.645186e-09 RADI step: 163 normalized residual: 8.264237e-09 RADI step: 164 normalized residual: 8.050445e-09 RADI step: 166 normalized residual: 6.138253e-09 RADI step: 167 normalized residual: 6.137569e-09 RADI step: 169 normalized residual: 5.414759e-09 RADI step: 170 normalized residual: 5.100608e-09 RADI step: 171 normalized residual: 5.098794e-09 RADI step: 173 normalized residual: 3.355112e-09 RADI step: 175 normalized residual: 3.178896e-09 RADI step: 177 normalized residual: 3.145711e-09 RADI step: 179 normalized residual: 4.807019e-10 RADI step: 181 normalized residual: 4.302671e-10 RADI step: 183 normalized residual: 4.127103e-10 RADI step: 185 normalized residual: 3.340906e-10 RADI step: 187 normalized residual: 3.001496e-10 RADI step: 189 normalized residual: 1.044920e-10 RADI step: 191 normalized residual: 9.084712e-11 mess_lrradi took 2.27 seconds outradi.Z: 3002 333

compare
if not(istest) figure(3); 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
