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, iSBN/ISSN: 978-3-901608-35-3. % % [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 program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, see <http://www.gnu.org/licenses/>. % % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % 2009-2020 %
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)=4.030239e+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'; tic; outnm = mess_lrnm(eqn, opts, oper); toc; if istest if min(outnm.res)>=opts.nm.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(outnm.res); title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('outnm.Z:'); disp(size(outnm.Z));
Using line search (res: 1.873798e+04) lambda: 3.444531e-03 NM step: 1 normalized residual: 9.731203e-01 relative change in K: 1.000000e+00 number of ADI steps: 221 Using line search (res: 2.675448e+01) lambda: 9.342475e-02 NM step: 2 normalized residual: 8.957650e-01 relative change in K: 1.455847e+00 number of ADI steps: 231 Using line search (res: 1.021836e+00) lambda: 5.879144e-01 NM step: 3 normalized residual: 2.623383e-01 relative change in K: 7.641741e-01 number of ADI steps: 260 NM step: 4 normalized residual: 2.290438e-02 relative change in K: 1.667444e-01 number of ADI steps: 275 NM step: 5 normalized residual: 3.166023e-04 relative change in K: 1.714730e-02 number of ADI steps: 244 NM step: 6 normalized residual: 5.116540e-07 relative change in K: 6.874317e-04 number of ADI steps: 257 NM step: 7 normalized residual: 1.532116e-12 relative change in K: 1.189840e-06 number of ADI steps: 240 Elapsed time is 13.725905 seconds. outnm.Z: 3002 332

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; tic; outradi = mess_lrradi(eqn, opts, oper); toc; if istest if min(outradi.res)>=opts.radi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(2); semilogy(outradi.res); title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM'); xlabel('number of iterations'); ylabel('normalized residual norm'); end disp('outradi.Z:'); disp(size(outradi.Z));
RADI step: 1 normalized residual: 4.729406e-01 RADI step: 2 normalized residual: 4.761866e-01 RADI step: 3 normalized residual: 5.698300e-01 RADI step: 4 normalized residual: 4.725082e-01 RADI step: 6 normalized residual: 4.711972e-01 RADI step: 8 normalized residual: 4.707539e-01 RADI step: 9 normalized residual: 4.641023e-01 RADI step: 11 normalized residual: 4.640827e-01 RADI step: 13 normalized residual: 4.640214e-01 RADI step: 14 normalized residual: 4.639910e-01 RADI step: 16 normalized residual: 4.636540e-01 RADI step: 18 normalized residual: 4.635961e-01 RADI step: 19 normalized residual: 4.506090e-01 RADI step: 21 normalized residual: 4.505054e-01 RADI step: 22 normalized residual: 4.504349e-01 RADI step: 24 normalized residual: 4.501128e-01 RADI step: 25 normalized residual: 4.500452e-01 RADI step: 26 normalized residual: 3.262328e-01 RADI step: 28 normalized residual: 3.256654e-01 RADI step: 30 normalized residual: 3.133978e-01 RADI step: 31 normalized residual: 3.133820e-01 RADI step: 32 normalized residual: 3.414981e-02 RADI step: 34 normalized residual: 3.396121e-02 RADI step: 36 normalized residual: 3.288124e-02 RADI step: 37 normalized residual: 3.288142e-02 RADI step: 38 normalized residual: 3.686013e-03 RADI step: 40 normalized residual: 3.304883e-03 RADI step: 42 normalized residual: 2.948404e-03 RADI step: 43 normalized residual: 2.864024e-03 RADI step: 44 normalized residual: 2.839482e-03 RADI step: 45 normalized residual: 2.839391e-03 RADI step: 47 normalized residual: 2.370904e-03 RADI step: 48 normalized residual: 2.325120e-03 RADI step: 49 normalized residual: 2.310415e-03 RADI step: 50 normalized residual: 2.310365e-03 RADI step: 52 normalized residual: 1.849386e-03 RADI step: 53 normalized residual: 1.832477e-03 RADI step: 55 normalized residual: 5.592773e-04 RADI step: 57 normalized residual: 5.592633e-04 RADI step: 58 normalized residual: 5.589192e-04 RADI step: 60 normalized residual: 3.026820e-04 RADI step: 62 normalized residual: 2.625747e-04 RADI step: 63 normalized residual: 2.527942e-04 RADI step: 64 normalized residual: 2.527736e-04 RADI step: 66 normalized residual: 2.064248e-04 RADI step: 68 normalized residual: 1.589102e-04 RADI step: 69 normalized residual: 1.523225e-04 RADI step: 70 normalized residual: 1.523196e-04 RADI step: 72 normalized residual: 1.452885e-04 RADI step: 74 normalized residual: 1.286848e-04 RADI step: 75 normalized residual: 1.276870e-04 RADI step: 77 normalized residual: 1.115358e-04 RADI step: 78 normalized residual: 1.115213e-04 RADI step: 80 normalized residual: 1.031334e-04 RADI step: 81 normalized residual: 1.018118e-04 RADI step: 82 normalized residual: 1.018071e-04 RADI step: 84 normalized residual: 3.896300e-05 RADI step: 86 normalized residual: 2.379412e-05 RADI step: 87 normalized residual: 2.039509e-05 RADI step: 89 normalized residual: 1.551376e-05 RADI step: 91 normalized residual: 1.063231e-05 RADI step: 92 normalized residual: 9.386292e-06 RADI step: 94 normalized residual: 8.836308e-06 RADI step: 96 normalized residual: 7.386663e-06 RADI step: 97 normalized residual: 6.769775e-06 RADI step: 99 normalized residual: 6.645817e-06 RADI step: 101 normalized residual: 5.064647e-06 RADI step: 103 normalized residual: 3.291714e-06 RADI step: 105 normalized residual: 1.210747e-06 RADI step: 107 normalized residual: 8.922208e-07 RADI step: 109 normalized residual: 4.698286e-07 RADI step: 111 normalized residual: 4.035859e-07 RADI step: 113 normalized residual: 3.574576e-07 RADI step: 115 normalized residual: 3.330885e-07 RADI step: 117 normalized residual: 2.759706e-07 RADI step: 119 normalized residual: 2.519761e-07 RADI step: 121 normalized residual: 2.445875e-07 RADI step: 123 normalized residual: 2.445763e-07 RADI step: 125 normalized residual: 2.408188e-07 RADI step: 127 normalized residual: 2.389848e-07 RADI step: 129 normalized residual: 2.373492e-07 RADI step: 131 normalized residual: 9.947402e-08 RADI step: 133 normalized residual: 2.549307e-08 RADI step: 135 normalized residual: 2.426164e-08 RADI step: 137 normalized residual: 1.326331e-08 RADI step: 139 normalized residual: 1.208237e-08 RADI step: 141 normalized residual: 1.184117e-08 RADI step: 143 normalized residual: 9.947308e-09 RADI step: 145 normalized residual: 9.496393e-09 RADI step: 147 normalized residual: 8.881583e-09 RADI step: 149 normalized residual: 5.975269e-09 RADI step: 151 normalized residual: 5.810398e-09 RADI step: 153 normalized residual: 4.893273e-09 RADI step: 155 normalized residual: 2.282152e-09 RADI step: 156 normalized residual: 2.201970e-09 RADI step: 158 normalized residual: 1.945668e-09 RADI step: 160 normalized residual: 1.431768e-09 RADI step: 161 normalized residual: 1.365336e-09 RADI step: 163 normalized residual: 1.070647e-09 RADI step: 165 normalized residual: 8.457081e-10 RADI step: 166 normalized residual: 8.119841e-10 RADI step: 168 normalized residual: 5.518412e-10 RADI step: 170 normalized residual: 4.458693e-10 RADI step: 171 normalized residual: 4.334268e-10 RADI step: 173 normalized residual: 3.378749e-10 RADI step: 175 normalized residual: 2.698316e-10 RADI step: 176 normalized residual: 2.640216e-10 RADI step: 178 normalized residual: 8.329984e-11 Elapsed time is 1.700862 seconds. outradi.Z: 3002 331

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-'); title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM'); xlabel('number of solves with A+p*M'); ylabel('normalized residual norm'); legend('LR-NM','RADI'); end
