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