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