Contents
function LQR_DAE3_SO(model,istest)
% Computes a Riccati feedback control for the constrained vibrational model % from [1] % % % Usage: % LQR_DAE3_SO(model,test) % % Input: % % model choice of the example system. Possible values % 'Stykel_small' Stykel's mass spring damper system from [1] % of dimension 600 in second order form (default) % 'Stykel_large' Stykel's mass spring damper system from [1] % of dimension 6000 in second order form % % 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] V. Mehrmann and T. Stykel. Balanced truncation model reduction for % large-scale systems in descriptor form. % In P. Benner, V. Mehrmann, and D. Sorensen, editors, Dimension % Reduction of Large-Scale Systems, volume 45 of Lecture Notes in % Computational Science and Engineering, pages 83–115. Springer-Verlag, % Berlin/Heidelberg, 2005. % 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,2); if (nargin < 1) model = 'Stykel_small'; end if ( nargin < 2 ) istest = 0; end
set operation
oper = operatormanager('dae_3_so');
load problem data
switch lower(model) case 'stykel_small' sys = load(sprintf('%s/../models/ms_ind3_by_t_stykel/g600.mat',... fileparts(mfilename('fullpath')))); case 'stykel_large' sys = load(sprintf('%s/../models/ms_ind3_by_t_stykel/g6000.mat',... fileparts(mfilename('fullpath')))); end eqn.M_=sys.M; eqn.E_=sys.D; eqn.K_=sys.K; eqn.G_=sys.G; eqn.haveE=1; eqn.alpha = -0.02; nv = size(eqn.M_,1); np = size(eqn.G_,1); eqn.B = sys.B(1:2*nv,:); eqn.C = sys.C(:,1:2*nv); eqn.type = 'T';
First we run the Newton-ADI Method
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); % 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 fprintf('2-Norm of the resulting feedback matrix: %g\n', norm(outnm.Z,2)); disp('size outnm.Z:'); disp(size(outnm.Z));
Warning: Non-stable Ritz values were detected. These will be removed from the set in further computations. NM step: 1 normalized residual: 1.415415e-01 relative change in K: 1.000000e+00 number of ADI steps: 46 Warning: Non-stable Ritz values were detected. These will be removed from the set in further computations. NM step: 2 normalized residual: 2.631219e-05 relative change in K: 1.382088e-02 number of ADI steps: 46 NM step: 3 normalized residual: 4.572588e-13 relative change in K: 1.822007e-06 number of ADI steps: 46 Elapsed time is 0.453428 seconds. 2-Norm of the resulting feedback matrix: 2.27051 size outnm.Z: 1200 50

Lets try the RADI metod and compare
RADI-MESS settings
opts.shifts.history = opts.shifts.num_desired*size(eqn.C,1); opts.shifts.num_desired=25; opts.shifts.num_Ritz=50; opts.shifts.num_hRitz=25; opts.shifts.b0=ones(2*nv+np,1); % 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.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: Non-stable Ritz values were detected. These will be removed from the set in further computations. RADI step: 1 normalized residual: 1.001499e+00 RADI step: 2 normalized residual: 4.884287e+00 RADI step: 4 normalized residual: 2.281476e+00 RADI step: 6 normalized residual: 3.593905e+00 RADI step: 8 normalized residual: 1.014958e+00 RADI step: 10 normalized residual: 1.425091e-01 RADI step: 12 normalized residual: 1.952333e-01 RADI step: 14 normalized residual: 7.656443e-02 RADI step: 16 normalized residual: 4.215812e-03 RADI step: 18 normalized residual: 2.983008e-03 RADI step: 20 normalized residual: 1.777674e-05 RADI step: 22 normalized residual: 7.529772e-07 RADI step: 24 normalized residual: 9.044695e-09 RADI step: 26 normalized residual: 5.459748e-10 RADI step: 27 normalized residual: 5.372400e-10 RADI step: 28 normalized residual: 9.524957e-10 RADI step: 30 normalized residual: 1.102023e-10 RADI step: 32 normalized residual: 1.669552e-10 RADI step: 34 normalized residual: 1.471539e-10 RADI step: 36 normalized residual: 1.558063e-11 Elapsed time is 0.191757 seconds. size outradi.Z: 1200 50

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
