Contents
function LQR_DAE2(problem,lvl,re,istest)
% Computes a stabilizing feedback by implicitly solving the generalized % Riccati equation for the equivalent projected system on the hidden manifold. % % Inputs: % problem either 'Stokes' or 'NSE' to choose the Stokes demo or the % linearized Navier-Stokes-Equation. (required) % % lvl discretization level 1 through 5 % (optional, only used in 'NSE' case, default: 1) % % re Reynolds number 300, 400, or 500 % (optional, only used in 'NSE' case, default: 500) % % istest flag to determine whether this demo runs as a CI test or % interactive demo % (optional, defaults to 0, i.e. interactive demo) % % Note that the 'NSE' option requires additional data available in a % separate 270MB archive and at least the 5th discretization level needs a % considerable amount of main memory installed in your machine. % % 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 %
Set operations
oper = operatormanager('dae_2');
Problem data
if nargin<1, problem='stokes'; end if nargin<2, lvl=1; end if nargin<3, re=500; end if nargin<4, istest=0; end switch lower(problem) case 'stokes' nin = 5; nout = 5; nx = 10; ny = 10; [eqn.E_,eqn.A_,eqn.B,eqn.C]=stokes_ind2(nin,nout,nx,ny); eqn.haveE=1; st=full(sum(diag(eqn.E_))); eqn.st=st; eqn.B=eqn.B(1:st,:); eqn.C=eqn.C(:,1:st); case 'nse' try load(sprintf('%s/../models/NSE/mat_nse_re_%d',... fileparts(mfilename('fullpath')),re),'mat'); catch error(['The files mat_nse_re_300.mat, mat_nse_re_400.mat and ', ... 'mat_nse_re_500.mat are available for dowload in a ', ... 'separate archive (270MB each). Please fetch them from the ', ... 'MESS download page and unpack them into the ', ... 'DEMOS/models/NSE folder.']); end eqn.A_=mat.mat_v.fullA{lvl}; eqn.E_=mat.mat_v.E{lvl}; eqn.haveE=1; eqn.B=mat.mat_v.B{lvl}; if re>200 opts.nm.K0=mat.mat_v.Feed_0{lvl}'; opts.radi.K0 = opts.nm.K0; end eqn.C=mat.mat_v.C{lvl}; eqn.st=mat.mat_mg.nv(lvl); otherwise error('input ''problem'' must be either ''NSE'' or ''Stokes'''); end
degrees of freedom: ------------------------------ total : 280 velocity : 180 pressure : 100 n_finite : 81 ------------------------------ Generating FVM matrices... -> Laplacians... -> gradient and divergence operator... -> B1, C1 and v1 velocity nodes... -> B2, C2 and v2 velocity nodes... Setting up system matrices ...
First we run the Newton-ADI Method
opts.norm = 2; % ADI tolerances and maximum iteration number opts.adi.maxiter = 300; opts.adi.res_tol = 1e-12; opts.adi.rel_diff_tol = 1e-16; opts.adi.info = 1; opts.adi.LDL_T=0; eqn.type='T';
n=size(eqn.A_, 1); opts.shifts.num_desired=5;%*nout; opts.shifts.num_Ritz=50; opts.shifts.num_hRitz=25; opts.shifts.method = 'projection'; opts.shifts.b0=ones(n,1);
Newton tolerances 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.projection.freq=0; opts.nm.projection.ortho=1; %opts.nm.projection.meth='care_nwt_fac'; opts.nm.res=struct('maxiter',10,'tol',1e-6,'orth',0); opts.nm.linesearch = 1; opts.nm.inexact = 'superlinear'; opts.nm.tau = 0.1; opts.nm.accumulateRes = 1;
use low-rank Newton-Kleinman-ADI
tic; outnm = mess_lrnm(eqn, opts, oper); toc; if not(istest) figure(1); disp(outnm.res); semilogy(outnm.res); title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM'); xlabel('number of newton iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outnm.Z:'); disp(size(outnm.Z));
ADI step: 1 normalized residual: 4.749604e-01 relative change in Z: 1.000000e+00 normalized outer residual: 4.748502e-01 ADI step: 2 normalized residual: 2.313480e-01 relative change in Z: 6.115557e-01 normalized outer residual: 2.309720e-01 ADI step: 3 normalized residual: 1.160895e-01 relative change in Z: 4.086903e-01 normalized outer residual: 1.154641e-01 ADI step: 4 normalized residual: 3.274932e-02 relative change in Z: 3.384842e-01 normalized outer residual: 3.186747e-02 ADI step: 5 normalized residual: 4.408328e-03 relative change in Z: 2.031583e-01 normalized outer residual: 3.477590e-03 outer tolerance: 1.188528e-02 NM step: 1 normalized residual: 3.477590e-03 relative change in K: 1.000000e+00 number of ADI steps: 5 ADI step: 1 normalized residual: 4.751369e-01 relative change in Z: 1.000000e+00 normalized outer residual: 4.746418e-01 ADI step: 2 normalized residual: 2.310230e-01 relative change in Z: 6.109172e-01 normalized outer residual: 2.308408e-01 ADI step: 3 normalized residual: 1.156490e-01 relative change in Z: 4.073025e-01 normalized outer residual: 1.155924e-01 ADI step: 4 normalized residual: 3.274385e-02 relative change in Z: 3.354002e-01 normalized outer residual: 3.273882e-02 ADI step: 5 normalized residual: 4.365684e-03 relative change in Z: 2.003457e-01 normalized outer residual: 4.365677e-03 ADI step: 6 normalized residual: 1.459721e-03 relative change in Z: 6.921311e-02 normalized outer residual: 1.459679e-03 ADI step: 7 normalized residual: 9.263121e-05 relative change in Z: 5.492722e-02 normalized outer residual: 9.253355e-05 outer tolerance: 3.863989e-04 NM step: 2 normalized residual: 9.253355e-05 relative change in K: 1.322216e-02 number of ADI steps: 7 ADI step: 1 normalized residual: 5.908896e-01 relative change in Z: 1.000000e+00 normalized outer residual: 5.902469e-01 ADI step: 2 normalized residual: 2.801793e-01 relative change in Z: 6.987959e-01 normalized outer residual: 2.799246e-01 ADI step: 3 normalized residual: 7.784030e-02 relative change in Z: 5.237500e-01 normalized outer residual: 7.780849e-02 ADI step: 4 normalized residual: 6.236314e-04 relative change in Z: 3.257720e-01 normalized outer residual: 6.234901e-04 ADI step: 5 normalized residual: 5.863780e-05 relative change in Z: 3.252883e-02 normalized outer residual: 5.863776e-05 ADI step: 6 normalized residual: 2.832583e-05 relative change in Z: 5.135531e-03 normalized outer residual: 2.832582e-05 ADI step: 7 normalized residual: 2.435343e-05 relative change in Z: 4.159017e-03 normalized outer residual: 2.435343e-05 ADI step: 8 normalized residual: 2.563350e-07 relative change in Z: 1.651384e-03 normalized outer residual: 2.563342e-07 outer tolerance: 3.304769e-06 NM step: 3 normalized residual: 2.563342e-07 relative change in K: 1.145646e-04 number of ADI steps: 8 ADI step: 1 normalized residual: 5.952306e-01 relative change in Z: 1.000000e+00 normalized outer residual: 5.945829e-01 ADI step: 2 normalized residual: 3.403191e-01 relative change in Z: 6.609782e-01 normalized outer residual: 3.399864e-01 ADI step: 3 normalized residual: 1.078645e-01 relative change in Z: 5.650374e-01 normalized outer residual: 1.078086e-01 ADI step: 4 normalized residual: 2.082073e-02 relative change in Z: 3.428347e-01 normalized outer residual: 2.081777e-02 ADI step: 5 normalized residual: 1.759708e-04 relative change in Z: 1.767522e-01 normalized outer residual: 1.759642e-04 ADI step: 6 normalized residual: 8.406516e-05 relative change in Z: 1.521109e-02 normalized outer residual: 8.406342e-05 ADI step: 7 normalized residual: 9.211191e-06 relative change in Z: 1.267083e-02 normalized outer residual: 9.211191e-06 ADI step: 8 normalized residual: 3.101855e-06 relative change in Z: 1.302219e-03 normalized outer residual: 3.101855e-06 ADI step: 9 normalized residual: 2.053261e-08 relative change in Z: 1.875062e-03 normalized outer residual: 2.053261e-08 ADI step: 10 normalized residual: 6.170187e-09 relative change in Z: 8.933795e-05 normalized outer residual: 6.170187e-09 ADI step: 11 normalized residual: 1.864742e-09 relative change in Z: 6.988242e-05 normalized outer residual: 1.864742e-09 ADI step: 12 normalized residual: 1.354434e-09 relative change in Z: 2.042627e-05 normalized outer residual: 1.354434e-09 ADI step: 13 normalized residual: 9.506925e-11 relative change in Z: 4.068725e-05 normalized outer residual: 9.506925e-11 outer tolerance: 3.943603e-09 NM step: 4 normalized residual: 9.506925e-11 relative change in K: 1.797500e-07 number of ADI steps: 13 Elapsed time is 0.115253 seconds. 0.0035 0.0001 0.0000 0.0000 size outnm.Z: 180 54

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 = opts.shifts.num_desired; % 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.shifts.info = 0; opts.radi.compute_sol_fac = 1; % Turned on for numerical stability reasons. opts.radi.get_ZZt = 0; opts.radi.maxiter = opts.adi.maxiter; 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 not(istest) figure(); semilogy(outradi.res); title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM'); xlabel('number of iterations'); ylabel('normalized residual norm'); end
RADI step: 1 normalized residual: 7.035353e-02 RADI step: 2 normalized residual: 1.197302e-02 RADI step: 3 normalized residual: 3.148166e-03 RADI step: 4 normalized residual: 2.370135e-03 RADI step: 5 normalized residual: 6.534372e-05 RADI step: 6 normalized residual: 5.589811e-05 RADI step: 7 normalized residual: 2.820919e-06 RADI step: 8 normalized residual: 2.015254e-07 RADI step: 9 normalized residual: 1.213884e-07 RADI step: 10 normalized residual: 7.433566e-09 RADI step: 11 normalized residual: 4.192365e-09 RADI step: 12 normalized residual: 7.152029e-10 RADI step: 13 normalized residual: 5.434915e-10 RADI step: 14 normalized residual: 4.511120e-10 RADI step: 15 normalized residual: 2.662895e-10 RADI step: 16 normalized residual: 2.161544e-10 RADI step: 17 normalized residual: 1.216147e-10 RADI step: 18 normalized residual: 6.488593e-11 Elapsed time is 0.083617 seconds.

compare
if istest if min(outnm.res)>=opts.nm.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end if min(outradi.res)>=opts.radi.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(); 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
