Поиск оптимального программного управления -- Sage

import time attach(DATA+"bvp.py") nf = lambda F: html('$$%s$$'%latex(F)) nfn = lambda F,N: html('$$%(n)s = %(f)s$$'%{'n':N, 'f':latex(F)}) nfp = lambda F,N,V: \ html(r'$$\frac{\partial %(n)s}{\partial %(v)s} = %(f)s$$'%{'n':N, 'f':latex(F), 'v':V})

aa = (9+1+9+1+0+7+9+5+3)/10+19/31 x1_0 = aa/3 x2_0 = -aa/5 t1 = 0 t2 = 10

var('x1 x2 u p1 p2 t')

(x1, x2, u, p1, p2, t)

(x1, x2, u, p1, p2, t)

f1 = x2 f2 = aa/8*x1 + u + sin(t)

H = -1/2*(3*x1^2 + 2*x2^2 + 3*u^2) + p1*f1 + p2*f2 nfn(H, 'H')

H = \frac{1}{1240} \, {\left(1240 \, u + 777 \, x_{1} + 1240 \, \sin\left(t\right)\right)} p_{2} + p_{1} x_{2} - \frac{3}{2} \, u^{2} - \frac{3}{2} \, x_{1}^{2} - x_{2}^{2}

H = \frac{1}{1240} \, {\left(1240 \, u + 777 \, x_{1} + 1240 \, \sin\left(t\right)\right)} p_{2} + p_{1} x_{2} - \frac{3}{2} \, u^{2} - \frac{3}{2} \, x_{1}^{2} - x_{2}^{2}

H.arguments()

(p1, p2, t, u, x1, x2)

(p1, p2, t, u, x1, x2)

eq_u = (H.differentiate(u)==0).simplify() nf(eq_u) res = solve(eq_u, u)[0] nf(res) u_at_max = res.rhs()

p_{2} - 3 \, u = 0
u = \frac{1}{3} \, p_{2}

p_{2} - 3 \, u = 0
u = \frac{1}{3} \, p_{2}

DHDP1 = H.differentiate(x1).simplify() DHDP2 = H.differentiate(x2).simplify() nfp(DHDP1, 'H', 'x_1') nfp(DHDP2, 'H', 'x_2')

\frac{\partial H}{\partial x_1} = \frac{777}{1240} \, p_{2} - 3 \, x_{1}
\frac{\partial H}{\partial x_2} = p_{1} - 2 \, x_{2}

\frac{\partial H}{\partial x_1} = \frac{777}{1240} \, p_{2} - 3 \, x_{1}
\frac{\partial H}{\partial x_2} = p_{1} - 2 \, x_{2}

eq1 = fast_callable(f1, vars=[x2], domain=RDF) eq2 = fast_callable(f2(u=u_at_max), vars=[x1, p2, t], domain=RDF) eq3 = fast_callable(DHDP1, vars=[x1, p2], domain=RDF) eq4 = fast_callable(DHDP2, vars=[x2, p1], domain=RDF)

func = lambda t, XP: \ [ eq1(XP[1]), \ eq2(XP[0],XP[3],t), \ eq3(XP[0],XP[3]), \ eq4(XP[1],XP[2]), \ ]

a = boundary_value_problem(func, [x1_0, x2_0, None, None], [None, None, 0, 0], t1, t2, [1,2,1,1], maxiter=1000) if a['error'] < 1.5: bvp_show_results(a, ['x_1', 'x_2', 'p_1', 'p_2', 'u_{max(H)}'], lambda Y: u_at_max(p2=Y[3]).n() )

Optimization terminated successfully.
         Current function value: 0.000042
         Iterations: 124
         Function evaluations: 224
Ошибка при вычислениях: 4.19146e-05
Начальный момент времени: 0.000000
Конечный момент времени: 10.000000
Значение функций в начальный момент времени:
 x_1   x_2   p_1   p_2   u_{max(H)}  
 259/155   -777/775   6.24749152923   -6.69962581765   -2.23320860588344  
Значения функций в конечный момент времени:
 x_1   x_2   p_1   p_2   u_{max(H)}  
 3.01170621575   3.28051313529   0.0023882783126   0.00601753822312   0.00200584607437199  
x_1: 

x_2: 

p_1: 

p_2: 

u_{max(H)}:

Optimization terminated successfully.
         Current function value: 0.000042
         Iterations: 124
         Function evaluations: 224
Ошибка при вычислениях: 4.19146e-05
Начальный момент времени: 0.000000
Конечный момент времени: 10.000000
Значение функций в начальный момент времени:
 x_1   x_2   p_1   p_2   u_{max(H)}  
 259/155   -777/775   6.24749152923   -6.69962581765   -2.23320860588344  
Значения функций в конечный момент времени:
 x_1   x_2   p_1   p_2   u_{max(H)}  
 3.01170621575   3.28051313529   0.0023882783126   0.00601753822312   0.00200584607437199  
x_1: 

x_2: 

p_1: 

p_2: 

u_{max(H)}:

Поиск оптимального программного управления

584 days ago by kpp