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