Scipy.optimize.minimize по оптимизации траектории (минимизация затрат) - PullRequest
1 голос
/ 10 марта 2019

Я отправил вопрос о том, как решить эту проблему в общем, здесь:

Оптимизация траектории для "Ракеты" с использованием scipy.optimize.minimize

В идеале я хотел бы просто минимизировать конечное время, но я не мог заставить оптимизатор добавлять время к переменной, которую можно корректно настроить, поэтому я решил, что вместо этого пока просто попытаюсь минимизировать u ^ 2. .

Вот код:


# Code
t_f = 1.0
t = np.linspace(0., t_f, num = 10) # Time array for 1 second into the future with 0.01 increment
u = np.zeros(t.size) + 650
print(u)
g = -650
initial_position = 0
initial_velocity = 0
final_position = 100
final_velocity = 100

def car_dynamics(x):
    # Create time vector
    # t = np.linspace(0., t_f, num = 100) # Time array for 1 second into the future with 0.01 increment


    # Integrate over entire time to find v as a function of t
    a = x + g
    v = int.cumtrapz(a, t, initial = 0) + initial_velocity

    # Integrate v(t) to get s(t)
    s = int.cumtrapz(v, t, initial = 0) + initial_position

    return s, v

def constraint1(x): # Final state constraints (Boundary conditions)
    s, v = car_dynamics(x)
    print('c1', s[0] - initial_position)
    return s[0] - initial_position

def constraint2(x): # Initial state constraints (initial conditions)
    s, v = car_dynamics(x)
    print('c2', v[0] - initial_velocity)
    return v[0] - initial_velocity

def constraint3(x):
    s, v = car_dynamics(x)
    print('c3', s[-1] - final_position)
    return s[-1] - final_position

def constraint4(x):
    s, v = car_dynamics(x)
    print('c4', v[-1] - final_velocity)
    return v[-1] - final_velocity

def constraint5(x):
    return x - 1000

def objective(x):
    u2 = np.square(x)
    return np.sum(u2)

cons = [{'type':'eq', 'fun':constraint1},
                {'type':'eq', 'fun':constraint2},
                {'type':'eq', 'fun':constraint3},
                {'type':'eq', 'fun':constraint4}]
                # {'type':'ineq', 'fun':constraint5}]


result = minimize(objective, u, constraints = cons, method = 'SLSQP', options={'eps':500, 'maxiter':1000, 'ftol':0.001, 'disp':True})
print(result)

Код выполняется, но оптимизатор не работает. Вот ошибка с выхода.

        c1 0.0
c2 0.0
c3 -100.0
c4 -100.0
c1 0.0
c2 0.0
c3 -100.0
c4 -100.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c3 -100.0
c3 -73.76543209876543
c3 -50.617283950617285
c3 -56.79012345679013
c3 -62.96296296296296
c3 -69.1358024691358
c3 -75.30864197530863
c3 -81.4814814814815
c3 -87.65432098765432
c3 -93.82716049382715
c3 -98.45679012345678
c4 -100.0
c4 -72.22222222222223
c4 -44.44444444444445
c4 -44.44444444444445
c4 -44.44444444444445
c4 -44.444444444444436
c4 -44.44444444444445
c4 -44.44444444444448
c4 -44.44444444444445
c4 -44.44444444444442
c4 -72.22222222222221
Singular matrix C in LSQ subproblem    (Exit mode 6)
            Current function value: 4225000.0
            Iterations: 1
            Function evaluations: 12
            Gradient evaluations: 1
     fun: 4225000.0
     jac: array([1800., 1800., 1800., 1800., 1800., 1800., 1800., 1800., 1800.,
       1800.])
 message: 'Singular matrix C in LSQ subproblem'
    nfev: 12
     nit: 1
    njev: 1
  status: 6
 success: False
       x: array([650., 650., 650., 650., 650., 650., 650., 650., 650., 650.])

Кажется, ограничения не соблюдаются в некотором количестве итераций. Должен ли я переключить свою целевую функцию, чтобы она содержала конечную скорость и конечную позицию? Я пробовал разные размеры шагов и что не с тем же кодом выхода.

Есть ли лучший способ использовать эту функцию для того, что я пытаюсь получить? Я пытаюсь получить управляющий вектор u (t) во всем интервале от t0 до t_f, чтобы я мог отправить эти команды ракете для оптимального управления. Сейчас я упростил оптимизацию до одной оси, просто чтобы узнать, как использовать функцию. Но, как вы видите, мне это не удалось.

Подобные примеры были бы чрезвычайно полезны, и я открыт для других методов оптимизации, если они являются числовыми и относительно быстрыми, поскольку я планирую в конечном итоге реализовать их в качестве прогнозирующего контроллера модели в реальном времени.

1 Ответ

2 голосов
/ 11 марта 2019

Ваша модель имеет как алгебраические, так и дифференциальные уравнения. Вам нужен решатель DAE для решения вышеуказанных неявных функций ODE. Один такой пакет, который мне известен, это гекко. (https://github.com/BYU-PRISM/GEKKO) Гекко специализируется на динамической оптимизации линейных, смешанных целочисленных и нелинейных задач оптимизации.

Ниже приведен пример проблемы запуска ракеты, которая сводит к минимуму время окончания. Доступный в http://apmonitor.com/wiki/index.php/Apps/RocketLaunch

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# create GEKKO model
m = GEKKO()

# scale 0-1 time with tf
m.time = np.linspace(0,1,101)

# options
m.options.NODES = 6
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.MAX_ITER = 500
m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0

# final time
tf = m.FV(value=1.0,lb=0.1,ub=100)
tf.STATUS = 1

# force
u = m.MV(value=0,lb=-1.1,ub=1.1)
u.STATUS = 1
u.DCOST = 1e-5

# variables
s = m.Var(value=0)
v = m.Var(value=0,lb=0,ub=1.7)
mass = m.Var(value=1,lb=0.2)

# differential equations scaled by tf
m.Equation(s.dt()==tf*v)
m.Equation(mass*v.dt()==tf*(u-0.2*v**2))
m.Equation(mass.dt()==tf*(-0.01*u**2))

# specify endpoint conditions
m.fix(s, pos=len(m.time)-1,val=10.0)
m.fix(v, pos=len(m.time)-1,val=0.0)

# minimize final time
m.Obj(tf)

# Optimize launch
m.solve()

print('Optimal Solution (final time): ' + str(tf.value[0]))

# scaled time
ts = m.time * tf.value[0]

# plot results
plt.figure(1)
plt.subplot(4,1,1)
plt.plot(ts,s.value,'r-',linewidth=2)
plt.ylabel('Position')
plt.legend(['s (Position)'])

plt.subplot(4,1,2)
plt.plot(ts,v.value,'b-',linewidth=2)
plt.ylabel('Velocity')
plt.legend(['v (Velocity)'])

plt.subplot(4,1,3)
plt.plot(ts,mass.value,'k-',linewidth=2)
plt.ylabel('Mass')
plt.legend(['m (Mass)'])

plt.subplot(4,1,4)
plt.plot(ts,u.value,'g-',linewidth=2)
plt.ylabel('Force')
plt.legend(['u (Force)'])

plt.xlabel('Time')
plt.show()

enter image description here

...