Contents

function LQR_TripleChain(n1, usfs, shifts, istest)
% LQR_TripleChain computes the optimal feedback for the LQR problem (e.g.
% [1]) with respect to the Triple-Chain-Oscillator from [2]. The
% computations use either first or second companion form linearization. In
% order to reduce the complexity, the 'so_1' and 'so_2' function handles
% cast all operations in the linearized 2n x 2n model back to operations
% with respect to the original n x n matrices [3,4]. The default function
% handles explicitly form the 2n x 2n matrices and treat them as a first
% order system.
%
% Usage:      LQR_TripleChain(n1, oper, shifts, istest)
%
% Inputs:
%
% n1          length of a single chain in the model
%             (optional; defaults to 1000)
%
% usfs        the set of user supplied functions to use.
%             possible values: 'so_1', 'so_2', 'default'
%             (optional; defaults to 'so_1'
%
% shifts      the desired ADI shift selection strategy
%             possible values: 'heur', 'projection'
%             (optional; defaults to 'projection'
%
% istest      flag to determine whether this demo runs as a CI test or
%             interactive demo
%             (optional, defaults to 0, i.e. interactive demo)
%
% References:
%
% [1] A. Locatelli, Optimal Control: An Introduction, Birkhäuser, Basel,
%     Switzerland, 2001.
%
% [2] 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
%
% [3] 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
%
% [4] P. Benner, P. Kürschner, J. Saak, An improved numerical method for
%     balanced truncation for symmetric second order systems,
%     Mathematical and Computer Modelling of Dynamical Systems 19 (6)
%     (2013) 593–615.
%     https://doi.org/10.1080/13873954.2013.794363

%
% 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)
%
narginchk(0, 4);
if nargin < 1
    n1 = 1000;
end
if nargin < 2
    usfs = 'so_1';
end
if nargin < 3
    shifts = 'projection';
end
if nargin < 4
    istest = false;
end

set operation

opts = struct();
[oper, opts] = operatormanager(opts, usfs);
% Problem data

alpha = 2;
Beta  = 5;
v     = 5;

switch usfs
    case 'so_1'
        [eqn.M_, eqn.E_, eqn.K_] = ...
            triplechain_MSD(n1, alpha, Beta, v);
        s  = size(eqn.K_, 1);
        eqn.B = [zeros(s, 1); ones(size(eqn.K_, 1), 1)];
        eqn.C = [ones(1, size(eqn.K_, 1)) zeros(1, s)];
    case 'so_2'
        [eqn.M_, eqn.E_, eqn.K_] = ...
            triplechain_MSD(n1, alpha, Beta, v);
        s  = size(eqn.K_, 1);
        eqn.B = [ones(size(eqn.K_, 1), 1); zeros(s, 1)];
        eqn.C = eqn.B';
    case 'default'
        [M_, E_, K_] = triplechain_MSD(n1, alpha, Beta, v);
        s = size(K_, 1);
        eqn.A_ = [sparse(s, s), -K_; -K_, -E_];
        eqn.E_ = [-K_, sparse(s, s); sparse(s, s), M_];
        eqn.B = [zeros(s, 1); ones(size(K_, 1), 1)];
        eqn.C = [ones(1, size(K_, 1)) zeros(1, s)];
        clear M_ E_ K_ s;
end

eqn.haveE = true;

ADI tolerances and maximum iteration number

opts.adi.maxiter          = 200;
opts.adi.res_tol          = 1e-10;
opts.adi.rel_diff_tol     = 0;
opts.adi.info             = 0;
opts.adi.accumulateK      = true;
opts.adi.accumulateDeltaK = false;
opts.adi.compute_sol_fac  = true;

eqn.type = 'T';

Heuristic shift parameters via basic Arnoldi

n = oper.size(eqn, opts);
switch shifts
    case 'heur'
        opts.shifts.num_desired = 25;
        opts.shifts.num_Ritz    = 50;
        opts.shifts.num_hRitz   = 25;

        opts.shifts.info        = 0;
        opts.shifts.method      = 'heur';
        opts.shifts.b0          = ones(n, 1);
    case 'projection'
        opts.shifts.num_desired = 6;
        opts.shifts.method      = 'projection';
        n                       = oper.size(eqn, opts);
        opts.shifts.b0          = ones(n, 1);
end
opts.shifts.truncate = 1e6; % remove all shifts larger than 1e6 or smaller
% than 1e-6 in absolute value in order to avoid
% loosing information about M or K in the
% shifted coefficients (p^2*M-pD+K)

Newton tolerances and maximum iteration number

opts.nm.maxiter                 = 25;
opts.nm.res_tol                 = 1e-10;
opts.nm.rel_diff_tol            = 1e-16;
opts.nm.info                    = 1;
opts.nm.accumulateRes           = true;
opts.nm.linesearch              = true;
opts.nm.inexact                 = 'quadratic';
opts.nm.tau                     = 0.1;
opts.norm                       = 'fro';
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(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
[mZ, nZ] = size(outnm.Z);
mess_fprintf(opts, 'size outnm.Z: %d x %d\n\n', mZ, nZ);
Warning: Could not compute initial projection shifts. Going to retry with random
right hand side. 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/TripleChain/LQR_TripleChain.m',156)">LQR_TripleChain (line 156)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lrnm.m',804)">mess_lrnm (line 804)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_para.m',379)">mess_para (line 379)</a>
  
		 Using line search (res: 2.241815e+04)  
		 lambda: 6.555361e-03  
  
 NM step:    1 normalized residual: 	2.694277e-01  
               relative change in K: 	1.000000e+00  
               number of ADI steps: 	2   
  
Warning: Could not compute initial projection shifts. Going to retry with random
right hand side. 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/TripleChain/LQR_TripleChain.m',156)">LQR_TripleChain (line 156)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lrnm.m',804)">mess_lrnm (line 804)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_para.m',379)">mess_para (line 379)</a>
  
 NM step:    2 normalized residual: 	1.000811e-02  
               relative change in K: 	9.991173e-02  
               number of ADI steps: 	3   
  
  
 NM step:    3 normalized residual: 	1.923044e-03  
               relative change in K: 	4.378956e-02  
               number of ADI steps: 	4   
  
  
 NM step:    4 normalized residual: 	3.062144e-04  
               relative change in K: 	1.748755e-02  
               number of ADI steps: 	5   
  
  
 NM step:    5 normalized residual: 	1.674489e-05  
               relative change in K: 	4.089949e-03  
               number of ADI steps: 	6   
  
  
 NM step:    6 normalized residual: 	6.314873e-08  
               relative change in K: 	2.511670e-04  
               number of ADI steps: 	26   
  
  
 NM step:    7 normalized residual: 	5.400188e-11  
               relative change in K: 	9.544375e-07  
               number of ADI steps: 	16   
  
mess_lrnm took   0.37 seconds   
size outnm.Z: 6002 x 1  
  

Lets try the RADI method and compare

RADI-MESS settings

opts.shifts.history = opts.shifts.num_desired * size(eqn.C, 1);

opts.shifts.method = 'projection';

% .. Suggest false (smart update is faster; convergence is the same).
opts.shifts.naive_update_mode  = false;
opts.radi.compute_sol_fac      = true;
opts.radi.get_ZZt              = true;
opts.radi.maxiter              = opts.adi.maxiter;
opts.norm                      = 2;
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(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
[mZ, nZ] = size(outradi.Z);
mess_fprintf(opts, 'size outradi.Z: %d x %d\n\n', mZ, nZ);
Warning: Could not compute initial projection shifts. Going to retry with random
right hand side. 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/TripleChain/LQR_TripleChain.m',191)">LQR_TripleChain (line 191)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lrradi.m',874)">mess_lrradi (line 874)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/shifts/mess_lrradi_get_shifts.m',105)">mess_lrradi_get_shifts (line 105)</a>
RADI step:    1 normalized residual: 9.906025e-01  
RADI step:    2 normalized residual: 5.645649e-01  
RADI step:    4 normalized residual: 8.970497e-03  
RADI step:    6 normalized residual: 4.171784e-07  
RADI step:    8 normalized residual: 3.505352e-10  
RADI step:   10 normalized residual: 2.016489e-10  
RADI step:   11 normalized residual: 1.179429e-10  
RADI step:   13 normalized residual: 1.014948e-10  
RADI step:   14 normalized residual: 9.926104e-11  
mess_lrradi took   0.12 seconds   
size outradi.Z: 6002 x 1  
  

compare

if istest
    nrm = norm(outnm.K - outradi.K, 'fro');
    nrmNM = norm(outnm.K, 'fro');
    if nrm / nrmNM >= 5e-9
        mess_err(opts, 'TEST:accuracy', ...
                 'unexpectedly inaccurate result: ', ...
                 '||K_NM - K_RADI||_F / ||K_NM||_F=%g', nrm / nrmNM);
    end
else
    figure(3);
    ls_nm   = [outnm.adi.niter];
    ls_radi = 1:outradi.niter;

    semilogy(cumsum(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