Tartalomjegyzék
\[ \newenvironment{dcases}{\left\{\begin{array}{ll}}{\end{array}\right.} \]Script d2018_03_18_without_centering_lv2_DT
Teljes Matlab script kiegészítő függvényekkel.
File: d2018_03_18_without_centering_lv2_DT.m Directory: /1_projects/2_sta/demonstrative Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. March 18.
Automatically generated stuff
Output:
┌d2018_03_18_without_centering_lv2_DT │ - Persistence for `d2018_03_18_without_centering_lv2_DT` reused (inherited) [run ID: 0509, 2018.03.18. Sunday, 15:21:05] │ - Script `d2018_03_18_without_centering_lv2_DT` backuped
U = @(A) A' + A;
n = 2;
P_generate_symvars_v1_caller(n,0);
A_lv = [
-2 -3
1.4 1
];
xeq = [1;1];
r = -A_lv*xeq;
fc_sym = diag(A_lv*x + r) * x;
% discretize:
h = 0.1;
f_sym = x + h*fc_sym;
f_fh = matlabFunction(f_sym, 'var', {t x});
Output:
│ ┌P_generate_symvars_v1_caller │ └ 0.23515 [sec]
Simulate
figure(6), hold on, axis([0 2 0 2]), axis equal, grid on;
plot(xeq(1),xeq(2),'o','MarkerSize',10,'LineWidth',4)
x0 = randn(n,1)*0.25 + xeq;
for i = 1:100
xuj = f_fh(0,x0);
% quiver(x0(1),x0(2),xuj(1)-x0(1),xuj(2)-x0(2));
plot([x0(1) xuj(1)],[x0(2) xuj(2)], 'k');
x0 = xuj;
end
Model representation
f0 = subs(f_sym, x, x*0);
f0_num = double(f0);
A_sym = P_kiemel(f_sym-f0,x);
pcz_symeq(A_sym*x + f0, f_sym, 'Kiemeles jol tortent: A(x)*x + f0 = f(x)');
warning off
lfr = P_LFR_reduction_v6(P_lfrdata_v6(sym2lfr(A_sym)),x);
warning on
[A,B,Pi] = deal(lfr.A,lfr.B,lfr.Pi);
Pib = [
1
x
Pi
];
Nb = P_annihilator_linsolve_v6(Pib,xd);
Pia = collect([
Pib
subs(Pi,x,f_sym)
],xd);
Pia_A_sym = P_kiemel(Pia(2:end),x);
warning off
Pia_lfr = P_LFR_reduction_v6(sym2lfr(Pia_A_sym),x);
warning on
H = blkdiag(1, [ Pia_lfr.A Pia_lfr.B ]);
tPia = [ 1 ; Pia_lfr.Pib ];
% H = eye(numel(Pia));
% tPia = Pia;
pcz_symeq(H*tPia,Pia, 'H*tPia = Pia')
Na = P_annihilator_linsolve_v6(tPia,xdDd);
p = numel(Pi);
[q,m] = size(Nb);
[qa,ma] = size(Na);
Ea = [ eye(n+p+1) zeros(n+p+1,p) ];
Aa = [ zeros(1,1+n+2*p) ; blkdiag([f0_num A B], eye(p)) ];
Aa(1) = 1;
P = sdpvar(m);
L = sdpvar(m,q);
La = sdpvar(ma,qa);
R = H'*(Aa'*P*Aa - Ea'*P*Ea)*H;
xlim = [
0.5 1.5
0.5 1.5
];
X_v = allcomb(xlim(1,:),xlim(2,:));
Nb_fh = matlabFunction(Nb,'vars',{x});
Na_fh = matlabFunction(Na,'vars',{x});
Optimization
CONS = [];
isposdef = @(P) P >= 0;
% isposdef = @(P) P >= eye(size(P))*0.0001;
for i = 1:size(X_v,1)
x_num = X_v(i,:)';
CONS = [
CONS
isposdef( P + L*Nb_fh(x_num) + Nb_fh(x_num)'*L' )
isposdef( -R + La*Na_fh(x_num) + Na_fh(x_num)'*La' )
];
end
sol = optimize(CONS,[],sdpsettings('solver','mosek','verbose',false));
pcz_info(sol.problem == 0, sol.info);
Visualization
P = value(P);
L = value(L);
R = value(R);
La = value(La);
tol = 1e-5;
isposdef = @(P) pcz_info(all(eig(P) > -tol), 'Positive semidefinite. Tolerance: %g',tol);
for i = 1:size(X_v,1)
x_num = X_v(i,:)';
isposdef( P + L*Nb_fh(x_num) + Nb_fh(x_num)'*L' )
% eig(P + L*Nb_fh(x_num) + Nb_fh(x_num)'*L')
isposdef( -R + La*Na_fh(x_num) + Na_fh(x_num)'*La' )
% eig(-R + La*Na_fh(x_num) + Na_fh(x_num)'*La')
end
resolution = 21;
[x1_num,x2_num] = meshgrid(linspace(xlim(1,1),xlim(1,2),resolution), linspace(xlim(2,1),xlim(2,2),resolution));
V = Pib'*P*Pib;
dV = collect(tPia'*R*tPia,x);
dV_elm = collect(subs(V,x,f_sym) - V,x);
pcz_symeq(dV_elm,dV,'Egyenlo-e a ket kifejezes dV-re?')
V_fh = matlabFunction(V);
dV_fh = matlabFunction(dV);
V_num = V_fh(x1_num,x2_num);
dV_num = dV_fh(x1_num,x2_num);
figure('Position', [ 572 305 1486 514 ], 'Color', [1 1 1]),
subplot(121), S = surf(x1_num,x2_num,V_num); axis vis3d tight; ptitle('$V(x)$')
S.CData = (S.CData).^0.025; colormap jet;
subplot(122), S = surf(x1_num,x2_num,dV_num); axis vis3d tight; ptitle('$\\Delta V(x) = V(x^+) - V(x)$')
S.CData = -(-S.CData).^0.25; colormap jet;
Output:
│ [ OK ] Kiemeles jol tortent: A(x)*x + f0 = f(x) │ [ OK ] The first n elements of sigma_b must be identity permutation (in theory) │ [ OK ] The upper left (n+k)x(n+k) submatrix of Theta should be full-rank │ [ OK ] Equation of the permuted model should be the same │ [ OK ] Equation of the reduced model should be the same │ ┌P_annihilator_linsolve_v6 │ │ ┌P_annihilator_linsolve_v6 - core computation │ │ └ 0.077842 [sec] │ │ ┌P_annihilator_linsolve_v6 - quick check I. - obtained annihilator, prec: 10 │ │ │ - maximal difference: 4.640385e-17, row nr. 7 │ │ └ 0.055347 [sec] │ │ ┌P_annihilator_linsolve_v6 - beautify annihilator + !roundings! │ │ └ 0.004683 [sec] │ │ ┌P_annihilator_linsolve_v6 - quick check II. - the beautified annihilator, prec: 10 │ │ │ - maximal difference: 3.364029e-16, row nr. 1 │ │ └ 0.037655 [sec] │ └ 0.17968 [sec] │ [ OK ] The first n elements of sigma_b must be identity permutation (in theory) │ [ OK ] The upper left (n+k)x(n+k) submatrix of Theta should be full-rank │ [ OK ] Equation of the permuted model should be the same │ [ OK ] Equation of the reduced model should be the same │ [ OK ] H*tPia = Pia │ ┌P_annihilator_linsolve_v6 │ │ ┌P_annihilator_linsolve_v6 - core computation │ │ └ 0.26232 [sec] │ │ ┌P_annihilator_linsolve_v6 - quick check I. - obtained annihilator, prec: 10 │ │ │ - maximal difference: 5.551115e-17, row nr. 1 │ │ └ 0.12228 [sec] │ │ ┌P_annihilator_linsolve_v6 - beautify annihilator + !roundings! │ │ └ 0.079731 [sec] │ │ ┌P_annihilator_linsolve_v6 - quick check II. - the beautified annihilator, prec: 10 │ │ │ - maximal difference: 4.420425e-16, row nr. 2 │ │ └ 0.04934 [sec] │ └ 0.51819 [sec] │ [ OK ] Successfully solved (MOSEK) │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Positive semidefinite. Tolerance: 1e-05 │ [ OK ] Egyenlo-e a ket kifejezes dV-re?
End of the script.
Output:
└ 5.5399 [sec]