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 ot the original nxn 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, Math. Comput. % Model. Dyn. Syst. 19 (6) (2013) 593–615. % https://doi.org/10.1080/13873954.2013.794363. % % 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 %
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=0; end
set operation
oper = operatormanager(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=1;
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 = 1;
opts.adi.accumulateDeltaK = 0;
opts.adi.compute_sol_fac = 1;
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 = 1; opts.nm.linesearch = 1; opts.nm.inexact = 'quadratic'; opts.nm.tau = 0.1; opts.norm = 'fro';
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); disp(outnm.res); 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('size outnm.Z:'); disp(size(outnm.Z));
Warning: Could not compute initial projection shifts. Going to retry with random right hand side. Using line search (res: 7.977985e+05) lambda: 1.052528e-03 NM step: 1 normalized residual: 4.681763e-01 relative change in K: 1.000000e+00 number of ADI steps: 3 Warning: Could not compute initial projection shifts. Going to retry with random right hand side. NM step: 2 normalized residual: 4.089779e-02 relative change in K: 1.981745e-01 number of ADI steps: 5 NM step: 3 normalized residual: 6.144845e-03 relative change in K: 7.811569e-02 number of ADI steps: 5 NM step: 4 normalized residual: 1.288759e-03 relative change in K: 3.585878e-02 number of ADI steps: 4 NM step: 5 normalized residual: 1.741846e-04 relative change in K: 1.319012e-02 number of ADI steps: 5 NM step: 6 normalized residual: 5.983070e-06 relative change in K: 2.444787e-03 number of ADI steps: 11 NM step: 7 normalized residual: 8.141213e-09 relative change in K: 9.018302e-05 number of ADI steps: 27 NM step: 8 normalized residual: 5.385013e-11 relative change in K: 1.230505e-07 number of ADI steps: 16 Elapsed time is 0.305592 seconds. 4.681762908708740e-01 4.089778698538684e-02 6.144844978242235e-03 1.288758669477708e-03 1.741845973617536e-04 5.983069943241979e-06 8.141212511790495e-09 5.385012678785875e-11 size outnm.Z: 6002 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 = 1; opts.radi.get_ZZt = 1; 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; 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('size outradi.Z:'); disp(size(outradi.Z));
Warning: Could not compute initial projection shifts. Going to retry with random right hand side. RADI step: 1 normalized residual: 9.797474e-01 RADI step: 2 normalized residual: 2.730705e-01 RADI step: 4 normalized residual: 7.516026e-04 RADI step: 6 normalized residual: 5.044813e-09 RADI step: 8 normalized residual: 4.046791e-10 RADI step: 10 normalized residual: 2.574211e-10 RADI step: 11 normalized residual: 1.321373e-10 RADI step: 13 normalized residual: 1.146582e-10 RADI step: 14 normalized residual: 1.128199e-10 RADI step: 15 normalized residual: 1.106538e-10 RADI step: 16 normalized residual: 1.073200e-10 RADI step: 17 normalized residual: 9.466406e-11 Elapsed time is 0.090930 seconds. size outradi.Z: 6002 1

compare
if istest nrm = norm(outnm.K-outradi.K,'fro'); nrmNM=norm(outnm.K,'fro'); if nrm/nrmNM >= 1e-9 error('MESS: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-'); 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
