2 %CPDE2016 This
function reproduces all the results from the CPDE paper
4 % Andreas Schmidt, 2016
6 % Function
for building all the results in the paper
11 % Reset the random number generator to
get the same results
18 % We consider a (damped) wave equation:
19 diffA = my_2d_matrix(n0, '0', '0', '0', '0.1');
20 A = @(mu) [zeros(n) eye(n); diffA -mu(1)*eye(n)];
21 B0 = [zeros(n,1); my_2d_vector(n0, '(0.3 <= x) & (x <= 0.7) & (0.4 <= y) & (y <= 0.6)')];
22 B = @(mu) B0/norm(B0);
23 Bw0 = [[zeros(n,1); my_2d_vector(n0, '(sin(pi*x)*sin(pi*y))')] zeros(N,1)];
24 Bw = @(mu) 20*Bw0/norm(Bw0);
26 Dz = @(mu) 0.1*[0; 1];
27 C0 = [my_2d_vector(n0, '(0.2 <= x) & (x <= 0.8) & (0.3 <= y) & (y <= 0.4)')' zeros(1,n)];
28 C = @(mu) C0/norm(C0);
29 Cz0 = [C0; zeros(1,N)];
30 Cz = @(mu) mu(2)*[Cz0/norm(Cz0)];
31 Dw = @(mu) 0.05*[0 1];
33 % Set up the RBmatlab models for the two algebraic Riccati equations:
34 descrP.A_comp = { [zeros(n) eye(n); diffA, zeros(n)], ...
35 blkdiag(zeros(n), -eye(n)) };
36 descrP.A_coeff = @(m)[ 1, m.mu(1) ];
37 descrP.B_comp = { B0/norm(B0) };
38 descrP.B_coeff = @(m) 1;
39 descrP.C_comp = { Cz0/norm(Cz0) };
40 descrP.C_coeff = @(m) m.mu(2);
41 descrP.R_comp = { (0.1*[0; 1])
'*(0.1*[0; 1]) };
43 descrP.mu_ranges = {[0.1, 2], [1, 100]};
44 descrP.mu_names = {'mu_damping
', 'mu_Q
'};
46 model_P = ARE.DefaultModel(descrP);
47 model_data_P = gen_model_data(model_P);
50 descrQ.A_comp = {descrP.A_comp{1}', descrP.A_comp{2}
'};
51 descrQ.B_comp = {C0'/norm(C0)};
52 descrQ.B_coeff = descrP.C_coeff;
53 descrQ.C_comp = { 20*Bw0
'/norm(Bw0) };
54 descrQ.C_coeff = @(m) 1;
55 descrQ.R_comp = { (0.05*[0 1])*(0.05*[0 1])' };
57 model_data_Q = gen_model_data(model_Q);
59 % Define the training parameters:
61 mu_train = gen_mus(N_train);
63 % Make sure we have a valid folder
for the results:
64 [~,~,~] = mkdir(
'result');
65 [~,~,~] = mkdir(
'plots');
67 % For the simulation we still need some more data:
68 t = linspace(0,5,500)
'; % The time span, final time T = 5
69 u = sin(t); % Just some input function
70 w = randn(length(t), 2); % Disturbance
71 x0 = [my_2d_vector(n0, '(sin(pi*x)*sin(pi*y))^2
'); zeros(n,1)];
76 for sk = 1:length(input)
79 display('Calculating all results
')
80 paper({'basis-Q
', 'basis-P
', 'basis-xhat
', 'basis-plot
', 'approaching-z2-plot
'})
83 display('Plotting results
for the basis generation plot
')
85 case 'full-dimensional
'
86 display('Performing a full dimensional simulation
')
88 % Assemble the full system:
89 sys = get_ss(mu_test);
90 % Perform uncontrolled simulation:
91 [sim_uc,~,x_uc] = lsim(sys, [0*u, w], t, [x0]);
93 % For comparison controlled simulation:
94 closed_loop = get_full_c_sys(mu_test);
95 [sim_c,~,x_c] = lsim(closed_loop, [0*u, w], t, [x0; 0*x0]);
97 plot(t, [sim_uc(:,1) sim_c(:,1)])
99 legend('Uncontrolled
', 'Controlled
')
101 case 'full-solution-3d
'
102 % Simulate the system with the noisy input and without any control.
103 % Output a surface plot
104 sys = get_ss([1 1 1]);
105 [~,~,x] = lsim(sys, [0*u, w], t, x0);
108 % Use the first n elements in the vector (they represent the
111 surf(linspace(0,1,n0),linspace(0,1,n0),reshape(X,n0,n0));
112 axis([0 1 0 1 min(min(x)) max(max(x))])
116 case 'uncontrolled_output_plot'
117 sys = get_ss([1 10 1]);
118 sim1 = lsim(sys, [0*u, w], t, [x0]);
119 sys2 = get_full_c_sys([1 10 1]);
120 sim2 = lsim(sys2, [0*u, w], t, [x0; 0*x0]);
122 p11 = plot(t, [sim1(:,3)], 'LineWidth', 0.8, 'Color', [0.65 0.65 0.65]);
123 plot(t, [sim2(:,3)], 'LineWidth', 0.8, 'Color', [0.65 0.65 0.65]);
124 p1 = plot(t, [sim1(:,4)], 'Color', 'b', 'LineWidth', 1.2);
125 p2 = plot(t, [sim2(:,4)], 'Color', 'r', 'LineWidth', 1.2);
126 legend([p1, p2], 'Uncontrolled', 'Controlled')
129 matlab2tikz( 'plots/output-example.tex', 'height', '\fheight', 'width', '\fwidth', 'parseStrings', false, 'extraAxisOptions','ylabel shift={-10pt}
' )
132 dsxhat = diag(Sxhat);
137 p = semilogy(r, dsp(r), '-
', r, dsq(r), '--
', r, dsxhat(r), '-.
');
138 p(1).LineWidth = [1.5];
139 p(2).LineWidth = [1.5];
140 p(3).LineWidth = [1.5];
142 legend('$\bar P$
', '$\bar Q$
', '$\bar x$
')
144 ylim = get(gca, 'YLim
')
145 plot(rp, dsp(rp), 'kd
', 'MarkerFaceColor
', 'k
')
146 line([rp rp], ylim, 'Color
', 'k
')
148 plot(rxhat, dsxhat(rxhat), 'kd
', 'MarkerFaceColor
', 'k
')
149 line([rxhat rxhat], ylim, 'Color
', 'k
')
151 plot(rq, dsq(rq), 'kd
', 'MarkerFaceColor
', 'k
')
152 line([rq rq], ylim, 'Color
', 'k
')
154 matlab2tikz( 'plots/svd-decay.tex
', 'height
', '\fheight
', 'width
', '\fwidth
', 'parseStrings
', false )
156 case 'reduced_full_simulation
'
158 mu_test = gen_mus(1);
160 full_closed_loop = get_full_c_sys(mu_test);
161 reduced_closed_loop = get_reduced_c_sys(mu_test, rp, rq, rxhat);
162 [full_sim_c,~,full_x_c] = lsim(full_closed_loop, [0*u, w], t, [x0; 0*x0]);
163 [reduced_sim_c,~,x_c] = lsim(reduced_closed_loop, [0*u, w], t, [x0; Vxhat(:,1:rxhat)'*0*x0]);
164 max(real(eig(reduced_closed_loop.a)))
165 plot(t, [full_sim_c(:,1), reduced_sim_c(:,1)])
169 case 'approaching-z2-plot'
170 this_t = linspace(0,1);
171 this_w = randn(100, 2);
172 this_u = randn(100,1);
174 %rxhat = size(Vxhat,2);
175 full_closed_loop = get_full_c_sys(mu_test);
176 this_u = 0*exp(2*this_t').^2;
177 [full_sim_c,~,full_x_c] = lsim(full_closed_loop, [this_u, this_w], this_t, [x0; 0*x0]);
180 plot(this_t, full_sim_c(:,2), 'linewidth', 1.2, 'DisplayName', 'Full');
182 sizes = [5 10 size(VP,2)];
183 ls = {
':',
'--',
'-.'};
185 for r = 1:length(sizes)
186 reduced_closed_loop = get_reduced_c_sys(mu_test, sizes(r), size(VQ,2), rxhat);
187 reduced_sim_c = lsim(reduced_closed_loop, [this_u, this_w], this_t, [x0; Vxhat(:,1:rxhat)'*0*x0]);
188 plot(this_t, reduced_sim_c(:,2), 'LineStyle', ls{r},
'LineWidth', 1);
190 legend({
'Full',
'N_P=5',
'N_P=10',
'N_P=60'},
'Location',
'southeast')
192 ylabel('Control u(t)')
193 matlab2tikz( 'plots/approaching_z2.tex', 'height', '\fheight', 'width', '\fwidth' )
197 display('Performing some performance checks')
207 p = care(full(a),b,c'*c,d'*d);
208 tsum = tsum + toc(f);
210 display(['Average time for solution of one
ARE is ', num2str(tsum/lmax)])
216 % Simulate the dynamical system
217 tmp_obs = ss(A(mut) - L*C(mut) - B(mut)*K, L, C(mut), []);
218 this_t = linspace(0,10);
219 this_y = sin(pi*this_t);
221 [asd,~,~] = lsim(tmp_obs, this_y, this_t);
224 V = Vxhat(:, 1:rxhat);
225 tmp_obs = ss(V'*(A(mut) - L*C(mut) - B(mut)*K)*V, V'*L, C(mut)*V, []);
227 [asd,~,~] = lsim(tmp_obs, this_y, this_t);
232 mu_test = gen_mus(N_test);
235 full_sys{k} = get_full_c_sys(mu_test(k,:));
238 NP = [5 10 20 40 60];
240 data = zeros(length(NX),length(NP));
241 grcs = @(m,np,nq,nx) get_reduced_c_sys(m,np,nq,nx);
244 parfor j = 1:length(NX)
248 rsys = grcs(mu_test(k,:),NP(i),2,NX(j));
250 fs = norm(full_sys{k}(1:2, 2:3),2);
251 r(k) = abs(norm(full_sys{k}(1:2, 2:3),2) - norm(rsys(1:2, 2:3), 2))/fs;
259 save('result/error_surface.mat', 'data');
261 case 'error-surface-output'
262 load('result/error_surface.mat');
265 NP = [5 10 20 40 60];
268 input.tableColLabels = arrayfun(@(i) num2str(NP(i)), 1:length(NP),
'UniformOutput',
false);
269 input.tableRowLabels = arrayfun(@(i) num2str(NX(i)), 1:length(NX),
'UniformOutput',
false);
270 input.dataFormat = {
'%.2e'};
273 % =======================================================================
274 % Functions
for the basis generation
275 % =======================================================================
277 % Implementation of the LRFG algorithm (not
using rb matlab
for
278 % portability). However,
this should be moved to RBmatlab at a certain
280 if exist(
'result/basis_P.mat',
'file')
281 display('Loading P basis from file')
282 load('result/basis_P.mat')
285 dd = gen_detailed_data(model_P, model_data_P);
286 [VP,P_info] = lrfg( @P, @Phat );
287 save('result/basis_P.mat', 'VP', 'P_info');
290 if exist('result/basis_Q.mat','file')
291 display('Loading Q basis from file')
292 load('result/basis_Q.mat')
295 [VQ,Q_info] = lrfg( @Q, @Qhat );
296 save('result/basis_Q.mat', 'VQ', 'Q_info');
299 display('Clearing the files')
300 delete('result/basis_P.mat', 'result/basis_Q.mat', ...
301 'result/basis_xhat.mat')
304 display('Loading all data matrices')
305 load('result/basis_Q.mat')
306 load('result/basis_P.mat')
307 load('result/basis_xhat.mat')
312 eE = sum(s(1:k).^2)/sum(s.^2);
320 display(' ... done loading the basis')
323 display('Performing basis generation for the regulator')
324 if exist('result/basis_xhat.mat')
325 load('result/basis_xhat.mat')
331 % Assemble the observer:
332 % sys = get_ss(mu_train(j,:));
333 display(['Iteration ', num2str(j)])
338 tmp_obs = ss(A(mut) - L*C(mut) - B(mut)*K, L, C(mut), []);
339 this_t = linspace(0,2*pi,500);
341 [~,~,x] = lsim(tmp_obs, this_y, this_t);
342 snapshots = [snapshots x'];
344 [Vxhat,Sxhat,~] = svd(snapshots);
346 save('result/basis_xhat.mat', 'Vxhat', 'Sxhat');
349 end % End of the "for"-loop
351 %% ========================================================================
352 function syscl = get_full_c_sys(mu)
353 % Get the system with full controller for comparison
354 syscl = get_reduced_c_sys(mu, [], [], []);
357 function syscl = get_reduced_c_sys(mu, red_p, red_q, red_xhat)
358 % Get the reduced systems. If one of the sizes red_p,red_q or red_xhat is
359 % empty, the full matrix/equation will be used.
363 if ~isempty(red_p) && red_p > 0
364 [~,khat] = detailedPhat(mu, red_p);
372 if ~isempty(red_q) && red_q > 0
373 [~,lhat] = Qhat(mu, red_q);
380 % 3. Assemble the controller:
381 if ~isempty(red_xhat) && red_xhat > 0
382 Vx = Vxhat(:, 1:red_xhat);
387 % Assemble the controller:
388 controller = ss(Vx'*(A(mu) - lhat*C(mu) - B(mu)*khat)*Vx, Vx'*lhat, -khat*Vx, []);
390 % Now close the loop and return the system:
391 % 1 --> u (Input of system)
392 % 3 --> y (Noisy output of system)
394 syscl = feedback(sys, controller, [1], [3], 1);
396 function p = gen_mus(N)
397 % Generate random parameter vectors (p in [4, 41])
398 p = repmat([1 1], N, 1) + repmat([3, 19], N, 1).*rand(N,2);
400 function sys = get_ss(mu)
401 % Assemble the system and return it
402 sys = ss(full(A(mu)), [B(mu) Bw(mu)], [Cz(mu); C(mu); C(mu)], [[Dz(mu)] ...
403 zeros(2); 0 Dw(mu); 0 0 0], 'inputname', {
'u', ...
404 'w1',
'w2'},
'outputname', {
'z1',
'z2',
'y',
'ytrue',});
407 function [p,k] = P(mu)
408 % This
function produces the
"true" feedback gain and the
true
409 % solution to the control
ARE
410 p = care(A(mu),B(mu),Cz(mu)
'*Cz(mu),Dz(mu)'*Dz(mu));
411 k = B(mu)
'*inv(Dz(mu)'*Dz(mu))*p;
414 function [phat,khat,PN,res] = Phat(mu,N,v)
415 % Calculate the reduced solution to a control
ARE
417 % Use the global VP matrix or the local one:
425 Qr = V'*Cz(mu)'*Cz(mu)*V;
428 PN = care(Ar, Br, Qr, R);
430 khat = B(mu)'*inv(R)*phat;
433 nPhat = norm(phat,'fro');
434 normA = norm(A(mu),'fro');
435 d = norm(Cz(mu)'*Cz(mu),'fro');
437 res = norm( A(mu)'*phat + phat*A(mu) - phat*B(mu)*inv(R)*B(mu)'*phat + Cz(mu)'*Cz(mu), 'fro') / (2*nPhat*normA + d + nPhat^2*norm(inv(R),'fro')*norm(B(mu),'fro')^2);
441 function [q,l] = Q(mu)
442 % Calculate the observer gain (full dimensional)
443 q = care(A(mu)', C(mu)', Bw(mu)*Bw(mu)', Dw(mu)*Dw(mu)');
444 l = q*C(mu)'*inv(Dw(mu)*Dw(mu)');
446 function [qhat,lhat,QN,res] = Qhat(mu,N,v)
447 % Calculate the reduced solution to a control
ARE (reduced)
454 QN = care(V'*A(mu)'*V, V'*C(mu)', V'*Bw(mu)*Bw(mu)'*V, Dw(mu)*Dw(mu)');
456 lhat = qhat*C(mu)'*inv(Dw(mu)*Dw(mu)');
458 nqhat = norm(qhat,'fro');
459 normA = norm(A(mu),'fro');
460 d = norm(Bw(mu)*Bw(mu)','fro');
461 res = norm( A(mu)*qhat + qhat*A(mu)' - qhat*C(mu)'*inv(Dw(mu)*Dw(mu)')*C(mu)*qhat + Bw(mu)*Bw(mu)', 'fro')/(d + 2*normA*nqhat + nqhat^2*norm(inv(Dw(mu)*Dw(mu)'),'fro')*norm(C(mu)','fro')^2);
465 % ================= Implementation of the greedy procedure ================
466 function [V,info] = lrfg( f, fr )
467 % Initialize the basis the most dominant mode of the first vector:
468 [V,~,~] = svds( f(mu_train(1,:)), 1);
470 errs = recalc_errs();
475 while max(errs) > tol
477 err_decay = [err_decay max(errs)];
479 display([num2str(i), ' - ', num2str(max(errs))])
481 x = f(mu_train(i,:));
483 % Mimick the low rank factor:
484 % Calculate the low rank factor:
487 x(:, find(diag(sqrt(d)) <= tol)) = [];
489 [u,s,~] = svd(x - V*V'*x, 'econ');
493 eE = sum(s(1:k))/sum(s);
498 display([' adding ', num2str(k), ' elements to the basis'])
499 V = orthonormalize_gram_schmidt([V, u(:,1:k)]);
500 errs = recalc_errs();
504 function errs = recalc_errs()
505 % Nested function that calculates the "error":
506 errs = zeros(size(mu_train,1),1);
507 parfor i = 1:size(mu_train,1)
508 [x,~,~,res] = fr(mu_train(i,:),[],V);
513 info.err_decay = err_decay;
518 function v = my_2d_vector(n0,f_str)
520 % Generates a vector v which contains the values of a function f(x,y)
521 % on an equidistant grid in the interior of the unit square. The grid
522 % points are numbered consistently with those used in the function
525 % This function is just used as an easy way to generate test problems
526 % rather than to solve PDEs.
530 % v = fdm_2d_vector( n0, f_str)
534 % n0 number of inner grid points in each dimension;
535 % f_str
string describing the function f in the space variables 'x'
536 % and 'y', e.g., f_str = 'sin(x+2*y)+3'.
540 % v vector of order n = n0^2 containing the values of f(x,y).
543 % LYAPACK 1.0 (Thilo Penzl, May 1999)
545 % Input data not completely checked!
550 error('Wrong number of input parameters.');
571 function [A, name] = my_2d_matrix(n0,fx_str,fy_str,g_str,mu_str)
573 % Generates the stiffness matrix A for the finite difference
574 % discretization (equidistant grid) of the PDE
576 % mu* laplace(u) - fx du/dx - fy du/dy - g u = r.h.s. on Omega
580 % Omega = (0,1)x(0,1) (unit square).
582 % This function is just used as an easy way to generate test problems
583 % rather than to solve PDEs.
587 % [A, name] = fdm_2d_matrix( n0, fx_str, fy_str, g_str )
591 % n0 number of inner grid points in each dimension;
592 % fx_str
string describing the function fx in the space variables
593 % 'x' and 'y', e.g., fx_str = 'sin(x+2*y)+3';
594 % fy_str
string describing the function fy in the space variables
596 % g_str
string describing the function g in the space variables
601 % A n-x-n sparse stiffness matrix, where n = n0^2;
602 % name
string describing the problem.
605 % LYAPACK 1.0 (Thilo Penzl, May 1999)
607 % Input data not completely checked!
615 name = ['FDM-2D: fx=',fx_str,'; fy=',fy_str,'; g=',g_str];
646 ptr = ptr+1; % A(i,i-n)
649 S(ptr) = mu*t2-fyv*t3;
653 ptr = ptr+1; % A(i,i-1)
656 S(ptr) = mu*t2-fxv*t3;
659 ptr = ptr+1; % A(i,i)
665 ptr = ptr+1; % A(i,i+1)
668 S(ptr) = mu*t2+fxv*t3;
672 ptr = ptr+1; % A(i,i+n0)
675 S(ptr) = mu*t2+fyv*t3;
681 A = -sparse(I,J,S,n2,n2);
function varargout = filecache_function(funcptr, varargin)
function used for file-caching other function calls.
function cpde2016(input)
CPDE2016 This function reproduces all the results from the CPDE paper.
Implementation of the parametric algebraic Riccati equation.