Tartalomjegyzék
\[ \newenvironment{dcases}{\left\{\begin{array}{ll}}{\end{array}\right.} \]Van der Pol egyszerűsített $\pi_a \left(x,\delta ,\delta^{˙} \right)$
Teljes Matlab script
(és live script)
kiegészítő függvényekkel.
Tekintsd meg LiveEditor nézetben is!
File: vdp_v7_egyszerusitett_Pia.m
Author: Peter Polcz (ppolcz@gmail.com)
Created on 2017. November 26.
% Automatically generated stuff
global SCOPE_DEPTH VERBOSE
SCOPE_DEPTH = 0;
VERBOSE = 0;
P_generate_symvars_v1_caller(2,0);
% define f(x):
e = 1;
A_sym = [
0 , -1
x1*x2 + 1 , -1
];
f_sym = expand(A_sym * x)
G_lfr = sym2lfr(A_sym);
lfr = P_lfrdata_v6(G_lfr);
lfr = P_LFR_reduction_v6(lfr, x);
Output:
f_sym =
-x2
x2*x1^2 + x1 - x2
A = lfr.A
B = lfr.B
Pi = lfr.Pi
Pib = lfr.Pib
n = numel(x);
p = numel(Pi);
m = numel(Pib);
Output:
A =
0 -1
1 -1
B =
0 0
1 0
Pi =
x1^2*x2
x1*x2
Pib =
x1
x2
x1^2*x2
x1*x2
dPi = expand(jacobian(Pi,x) * f_sym)
Pia = [Pib ; dPi]
Pia_A_sym = P_kiemel(Pia,x)
Pia_lfr = P_LFR_reduction_v6(sym2lfr(Pia_A_sym),x);
H = [ Pia_lfr.A Pia_lfr.B ]
tPia = Pia_lfr.Pib
Output:
dPi =
x1^4*x2 + x1^3 - x1^2*x2 - 2*x1*x2^2
x1^3*x2 + x1^2 - x1*x2 - x2^2
Pia =
x1
x2
x1^2*x2
x1*x2
x1^4*x2 + x1^3 - x1^2*x2 - 2*x1*x2^2
x1^3*x2 + x1^2 - x1*x2 - x2^2
Pia_A_sym =
[ 1, 0]
[ 0, 1]
[ x1*x2, 0]
[ x2, 0]
[ x1^3*x2 + x1^2 - x1*x2 - 2*x2^2, 0]
[ x2*x1^2 + x1 - x2, -x2]
H =
1 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0
0 0 -1 1 0 1 0 0 2 0
0 0 0 0 1 0 1 -1 0 1
tPia =
x1
x2
x1^2*x2
x1^4*x2
x1^3*x2
x1^3
x1^2
x1*x2
-x1*x2^2
-x2^2
Nb_sym = P_annihilator_linsolve_v6(Pib,x)
Na_sym = P_annihilator_linsolve_v6(tPia,x)
Na = matlabFunction(Na_sym,'vars',{x});
Nb = matlabFunction(Nb_sym,'vars',{x});
Output:
Nb_sym = [ x2, 0, 0, -1] [ 0, x1, 0, -1] [ 0, 0, 1, -x1] Na_sym = [ x1, 0, 0, 0, 0, 0, -1, 0, 0, 0] [ x2, 0, 0, 0, 0, 0, 0, -1, 0, 0] [ 0, x1, 0, 0, 0, 0, 0, -1, 0, 0] [ 0, x2, 0, 0, 0, 0, 0, 0, 0, 1] [ 0, 0, 1, 0, 0, 0, 0, -x1, 0, 0] [ 0, 0, x1, 0, 0, -x2, 0, 0, 0, 0] [ 0, 0, x2, 0, 0, 0, 0, 0, x1, 0] [ 0, 0, 0, 1, -x1, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, 1, -x2, 0, 0, 0, 0] [ 0, 0, 0, 0, 0, 1, -x1, 0, 0, 0] [ 0, 0, 0, 0, 0, 0, x2, -x1, 0, 0] [ 0, 0, 0, 0, 0, 0, 0, x2, 0, x1] [ 0, 0, 0, 0, 0, 0, 0, 0, 1, -x1]
Aa = [
A B zeros(n,p)
zeros(p,m) eye(p)
];
P = sdpvar(m);
Pa = [ eye(m) ; zeros(p,m) ] * P * Aa;
R = H' * (Pa' + Pa) * H;
La = sdpvar(size(Na_sym,2),size(Na_sym,1),'full');
Lb = sdpvar(size(Nb_sym,2),size(Nb_sym,1),'full');
U = @(A) A + A';
X_v = [
-1.1463 -2.3221
-0.2645 -2.3221
1.7930 -0.2645
1.7930 0.9112
1.2051 2.3809
0.3233 2.3809
-1.7343 0.3233
-1.7343 -1.1463
-1.5227 -1.6840
-1.7950 -0.2723
-0.7462 -2.4000
-1.8555 -0.2218
0.8112 2.5000
-0.7296 -2.5204
-1.6480 -1.5306
-1.9133 -0.5306
1.9133 0.2449
1.5561 1.7245
1.6582 1.8469
-1.2194 0.8878
-1.3716 0.7544
1.3747 -0.7296
] * 0.8;
chull = convhull(X_v(:,1),X_v(:,2));
xv = X_v(chull,1);
yv = X_v(chull,2);
figure, plot(xv,yv), grid on
Constraints = [];
for i = 1:size(X_v,1)
w = X_v(i,:)';
Constraints = [ Constraints ...
P + U(Lb*Nb(w)) >= 0
R + U(La*Na(w)) <= 0
];
end
sol = optimize(Constraints)
Output:
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 152
Cones : 0
Scalar variables : 0
Matrix variables : 44
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator - tries : 0 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 152
Cones : 0
Scalar variables : 0
Matrix variables : 44
Integer variables : 0
Optimizer - threads : 4
Optimizer - solved problem : the primal
Optimizer - Constraints : 152
Optimizer - Cones : 0
Optimizer - Scalar variables : 0 conic : 0
Optimizer - Semi-definite variables: 44 scalarized : 1430
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 1.01e+04 after factor : 1.01e+04
Factor - dense dim. : 0 flops : 2.99e+06
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 2.2e+01 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.01
1 3.7e+00 1.7e-01 4.0e-01 1.00e+00 0.000000000e+00 0.000000000e+00 1.7e-01 0.04
2 7.0e-01 3.2e-02 1.8e-01 1.00e+00 0.000000000e+00 0.000000000e+00 3.2e-02 0.05
3 3.2e-02 1.4e-03 3.7e-02 1.00e+00 0.000000000e+00 0.000000000e+00 1.4e-03 0.06
4 2.0e-06 9.0e-08 2.9e-04 1.00e+00 0.000000000e+00 0.000000000e+00 9.0e-08 0.07
5 2.1e-12 9.6e-14 9.4e-14 1.00e+00 0.000000000e+00 0.000000000e+00 9.4e-14 0.08
Optimizer terminated. Time: 0.19
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 0.0000000000e+00 nrm: 2e-11 Viol. con: 6e-12 barvar: 0e+00
Dual. obj: 0.0000000000e+00 nrm: 2e+00 Viol. con: 0e+00 barvar: 1e-13
Optimizer summary
Optimizer - time: 0.19
Interior-point - iterations : 5 time: 0.08
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00
sol =
struct with fields:
yalmiptime: 0.3387
solvertime: 0.2033
info: 'Successfully solved (MOSEK)'
problem: 0
sol
check(Constraints)
P = value(P);
R = value(R);
V_sym = Pib'*P*Pib;
% dV_sym = tPia'*R*tPia;
dV_sym = jacobian(V_sym,x)*f_sym;
V = matlabFunction(V_sym);
dV = matlabFunction(dV_sym);
xmin = min(X_v(:,1));
xmax = max(X_v(:,1));
ymin = min(X_v(:,2));
ymax = max(X_v(:,2));
res = 100;
xres = round((xmax - xmin) * res);
yres = round((ymax - ymin) * res);
[xx1,xx2] = meshgrid(linspace(xmin,xmax,xres), linspace(ymin,ymax,yres));
V_num = V(xx1,xx2);
dV_num = dV(xx1,xx2);
figure('Position', [ 276 433 1331 416 ], 'Color', [1 1 1]),
subplot(121), hold on
surf(xx1,xx2,V_num), shading interp; light
contour(xx1,xx2,V_num)
view([ -37.5 , 30 ])
subplot(122),
surf(xx1,xx2,dV_num), shading interp; light
Output:
sol =
struct with fields:
yalmiptime: 0.3387
solvertime: 0.2033
info: 'Successfully solved (MOSEK)'
problem: 0
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Primal residual| Dual residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Matrix inequality| 0.44015| 1.9578e-14|
| #2| Matrix inequality| 0.072979| 1.6942e-14|
| #3| Matrix inequality| 0.55009| 1.6896e-14|
| #4| Matrix inequality| 0.13522| 1.3339e-14|
| #5| Matrix inequality| 0.332| 8.8256e-15|
| #6| Matrix inequality| 0.021842| 1.0614e-14|
| #7| Matrix inequality| 0.33555| 8.0214e-15|
| #8| Matrix inequality| 0.030098| 1.1553e-14|
| #9| Matrix inequality| 0.42532| 6.5821e-15|
| #10| Matrix inequality| 0.078285| 1.2779e-14|
| #11| Matrix inequality| 0.57217| 7.0473e-15|
| #12| Matrix inequality| 0.11565| 1.077e-14|
| #13| Matrix inequality| 0.39507| 1.3622e-14|
| #14| Matrix inequality| 0.034383| 9.7487e-15|
| #15| Matrix inequality| 0.38119| 1.638e-14|
| #16| Matrix inequality| 0.030096| 1.3266e-14|
| #17| Matrix inequality| 0.40875| 1.7985e-14|
| #18| Matrix inequality| 0.044486| 1.6769e-14|
| #19| Matrix inequality| 0.38051| 1.4766e-14|
| #20| Matrix inequality| 0.071962| 1.0628e-14|
| #21| Matrix inequality| 0.49314| 1.8673e-14|
| #22| Matrix inequality| 0.13396| 1.5077e-14|
| #23| Matrix inequality| 0.36393| 1.473e-14|
| #24| Matrix inequality| 0.054104| 1.0389e-14|
| #25| Matrix inequality| 0.49148| 6.5182e-15|
| #26| Matrix inequality| 0.12377| 1.2858e-14|
| #27| Matrix inequality| 0.48355| 1.894e-14|
| #28| Matrix inequality| 0.11603| 1.4469e-14|
| #29| Matrix inequality| 0.38771| 1.7573e-14|
| #30| Matrix inequality| 0.022603| 1.4861e-14|
| #31| Matrix inequality| 0.34639| 1.521e-14|
| #32| Matrix inequality| 0.020188| 1.1532e-14|
| #33| Matrix inequality| 0.30716| 8.2921e-15|
| #34| Matrix inequality| 0.024224| 1.0083e-14|
| #35| Matrix inequality| 0.37521| 7.3749e-15|
| #36| Matrix inequality| 0.047215| 1.2066e-14|
| #37| Matrix inequality| 0.34966| 7.1193e-15|
| #38| Matrix inequality| 0.012399| 5.6006e-23|
| #39| Matrix inequality| 0.52729| 1.1714e-14|
| #40| Matrix inequality| 0.065655| 1.0137e-14|
| #41| Matrix inequality| 0.48931| 1.2216e-14|
| #42| Matrix inequality| 0.046868| 9.252e-15|
| #43| Matrix inequality| 0.42941| 1.0083e-14|
| #44| Matrix inequality| 0.046219| 1.0911e-14|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
figure('Position', [ 276 433 1331 416 ], 'Color', [1 1 1]),
V_num(inpolygon(xx1,xx2,xv,yv) == 0) = NaN;
dV_num(inpolygon(xx1,xx2,xv,yv) == 0) = NaN;
subplot(121), hold on
surf(xx1,xx2,V_num), shading interp; light
contour(xx1,xx2,V_num)
view([ -37.5 , 30 ])
subplot(122),
surf(xx1,xx2,dV_num), shading interp; light
Eddig volt: 344
nr_Variables = ...
m*(m+1)/2 + ... P
numel(Lb) + numel(La)
Output:
nr_Variables = 152