Блог

Всех с праздником! У меня наконец-то появилось свободное время для работы над моим проектом, и, конечно, я застрял, как обычно, лол.

Я ищу рекомендации / примеры, которые позволили бы мне смоделировать следующее:

1. У меня есть входные данные (назовем их «переход»), которые являются двоичными (0 или 1), и я хочу, чтобы их можно было использовать только один раз (или, возможно, «n» количество раз, где n <# временных шагов) за весь временной горизонт. Влияние, которое этот ввод оказывает на систему, заключается в том, что он мгновенно увеличивает скорость системы на некоторую заранее определенную величину.
2. Второе влияние, которое это оказывает на систему, заключается в том, что изменяется набор динамик, которые продвигают систему вперед во времени. (В этом случае «скачок» изменит динамику от движущей системы к летающей системе). В будущем также будет «двойной скачок», который не меняет динамику, но все же обеспечивает мгновенное изменение скорости. В настоящее время я пытаюсь закончить первую часть, а затем попытаюсь реализовать это. Просто хочу, чтобы мое видение было понятным для всех, кто это читает.
3. Также еще одна часть, предназначенная для будущего модели: я хотел бы иметь возможность взаимодействовать с системой с шаровым объектом, скажем, с помощью if2 / if3, и если положение системы находится в некотором радиусе от положения другого объекта, импульс будет передаваться на зависящий от шарового объекта объекто таких вещах, как скорости шара и системы. Чтобы все было правильно, я полагаю, мне нужен способ определить временной шаг, который происходит в точке взаимодействия, что, я полагаю, означает, что мне понадобится какой-то переменный временной вектор. Любые примеры для них будут высоко оценены.

Хорошо, так что 2 и 3 здесь только для того, чтобы быть здесь, а не для того, чтобы быть основными моментами этого вопроса. Я думаю, что смогу разобраться с ними, как только смогу разобраться с реализацией этого странного ввода «скачка».

Мой текущий план состоит в том, чтобы иметь MV с именем ‘u_jump’, который не является целым числом. Затем есть переменная с именем ‘jump_hist’, которая по сути является «интегралом» от ‘u_jump’, и я даю jump_hist верхнюю границу 1. Что я делаю прямо сейчас, так это просто притворяюсь, что этот u_jump является ускорением в системе путем добавления к уравнению velocity.dt() . Это работает теоретически, но на самом деле не представляет систему, которую я пытаюсь контролировать идеально.

Какой был бы лучшим примером для меня, чтобы извлечь некоторые уроки для реализации этого? И еще один вопрос: есть ли способ заставить решатель IPOPT работать для целочисленных решений, задав целым числам допуск? Что-то вроде опции minlp solver ‘minlp_integer_tol 0.05’, таким образом, я все еще могу получить скорость IPOPT, но возможность включать переменные / уравнения целочисленного стиля, такие как if3 () и т.д… Если нет, есть ли способы, которыми я могу приблизиться к целочисленному решению с помощью нецелочисленного решения, чтобы при реализации управления в реальной системе разница между нецелочисленным решением и целочисленным решением находилась в пределах некоторого приемлемого допуска, чтобы считать это помехой, которую контроллер обратной связи мог бы смягчить?

Я знаю, что мои вопросы всегда хаха. Надеюсь, это будет полезно для других в будущем! Вот мой код на данный момент. Дайте мне знать, если код вызывает у вас проблемы или что-нибудь, что я мог бы прояснить в вопросе.

О, и последнее замечание, в настоящее время это настроено как 2D-система полета. Я удалил динамику движения (закомментированные сплайны c) для простоты реализации этого ввода «скачка».

Еще раз всех с праздниками!

 import numpy as np
import matplotlib.pyplot as plt
import math
import gekko
from gekko import GEKKO
import csv
from mpl_toolkits.mplot3d import Axes3D

class Optimizer():
    def __init__(self):
#################GROUND DRIVING OPTIMIZER SETTTINGS##############
        self.d = GEKKO(remote=False) # Driving on ground optimizer

        ntd = 21

        self.d.time = np.linspace(0, 1, ntd) # Time vector normalized 0-1

        # options
        self.d.options.NODES = 3
        self.d.options.SOLVER = 3
        self.d.options.IMODE = 6# MPC mode
        self.d.options.MAX_ITER = 800
        self.d.options.MV_TYPE = 0
        self.d.options.DIAGLEVEL = 0
        # self.d.options.OTOL = 1

        # final time for driving optimizer
        self.tf = self.d.FV(value=1.0,lb=0.1,ub=10.0, name='tf')

        # allow gekko to change the tf value
        self.tf.STATUS = 1
        
        # time variable
        self.t = self.d.Var(value=0)
        self.d.Equation(self.t.dt()/self.tf == 1)

        # Acceleration variable
        self.a = self.d.MV(fixed_initial=False, lb = 0, ub = 1, name='a')
        self.a.STATUS = 1
        
        # Jumping integer varaibles and equations
        self.u_jump = self.d.MV(fixed_initial=False, lb=0, ub=1, integer=True)
        self.u_jump.STATUS = 1
        self.jump_hist = self.d.Var(value=0, name='jump_hist', lb=0, ub=1)
        self.d.Equation(self.jump_hist.dt() == self.u_jump*(ntd-1))
        # self.d.Equation(1.0 >= self.jump_hist)

        # pitch input throttle (rotation of system)
        self.u_p = self.d.MV(fixed_initial=False, lb = -1, ub=1)
        self.u_p.STATUS = 1

        # Final variable that allows you to set an objective function considering only final state
        self.p_d = np.zeros(ntd)
        self.p_d[-1] = 1.0
        self.final = self.d.Param(value = self.p_d, name='final')

        # Model constants and parameters
        self.Dp = self.d.Const(value = 2.7982, name='D_pitch')
        self.Tp = self.d.Const(value = 12.146, name='T_pitch')
        self.pi = self.d.Const(value = 3.14159, name='pi')
        self.g = self.d.Const(value = 0, name='Fg')
        self.jump_magnitude = self.d.Param(value = 3000, name = 'jump_mag')

    def optimize2D(self, si, sf, vi, vf, ri, omegai): #these are 1x2 vectors s or v [x, z]

        # variables and intial conditions 
        # Position in 2d
        self.sx = self.d.Var(value=si[0], lb=-4096, ub=4096, name='x') #x position
        # self.sy = self.d.Var(value=si[1], lb=-5120, ub=5120, name='y') #y position
        self.sz = self.d.Var(value = si[1])

        # Pitch rotation and angular velocity
        self.pitch = self.d.Var(value = ri, name='pitch', lb=-1*self.pi, ub=self.pi)
        self.pitch_dot = self.d.Var(fixed_initial=False, name='pitch_dot')

        # Velocity in 2D
        self.v_mag = self.d.Var(value=(vi), name='v_mag')
        self.vx = self.d.Var(value=np.cos(ri), name='vx') #x velocity
        # self.vy = self.d.Var(value=(np.sin(ri) * vi), name='vy') #y velocity
        self.vz = self.d.Var(value = (np.sin(ri) * vi), name='vz')

## Non-linear state dependent dynamics descired as csplines.
        #curvature vs vel as a cubic spline for driving state
        cur = np.array([0.0069, 0.00398, 0.00235, 0.001375, 0.0011, 0.00088])
        v_cur = np.array([0,500,1000,1500,1750,2300])
        v_cur_fine = np.linspace(0,2300,100)
        cur_fine = np.interp(v_cur_fine, v_cur, cur)
        self.curvature = self.d.Var(name='curvature')
        self.d.cspline(self.v_mag, self.curvature, v_cur_fine, cur_fine)

        # throttle vs vel as cubic spline for driving state
        ba=991.666 #Boost acceleration magnitude
        kv = np.array([0, 1410, 2300]) #velocity input
        ka = np.array([1600 ba, 0 ba, 0 ba]) #acceleration ouput
        kv_fine = np.linspace(0, 2300, 100) # Higher resolution
        ka_fine = np.interp(kv_fine, kv, ka) # Piecewise linear high resolution of ka
        self.throttle_acceleration = self.d.Var(fixed_initial=False, name='throttle_accel')
        self.d.cspline(self.v_mag, self.throttle_acceleration, kv_fine, ka_fine)

# Differental equations
        # Velocity diff eqs
        self.d.Equation(self.vx.dt()/self.tf == (self.a*ba * self.d.cos(self.pitch)*self.jump_hist)   (self.a * self.throttle_acceleration * (1-self.jump_hist))   (self.u_jump * self.jump_magnitude * self.d.cos(self.pitch   np.pi/2)))
        self.d.Equation(self.vz.dt()/self.tf == (self.a*ba * self.d.sin(self.pitch)*self.jump_hist) - (self.g * (1-self.jump_hist))   (self.u_jump * self.jump_magnitude * self.d.sin(self.pitch   np.pi/2)))
        self.d.Equation(self.v_mag == self.d.sqrt((self.vx*self.vx)   (self.vz*self.vz)))
        self.d.Equation(2300 >= self.v_mag)

        # Position diff eqs
        self.d.Equation(self.sx.dt()/self.tf == self.vx)
        # self.d.Equation(self.sy.dt()/self.tf == self.vy)
        self.d.Equation(self.sz.dt()/self.tf == self.vz)

        # Orientation diff eqs
        self.d.Equation(self.pitch_dot.dt()/self.tf == ((self.Tp * self.u_p)   (self.Dp * self.pitch_dot * (1 - self.d.abs2(self.u_p)))) * self.jump_hist)
        self.d.Equation(self.pitch.dt()/self.tf == self.pitch_dot)


# Objective functions
        # Final Position Objectives
        self.d.Minimize(self.final*1e2*((self.sz-sf[1])**2)) # z final position objective
        self.d.Minimize(self.final*1e2*((self.sx-sf[0])**2)) # x final position objective
        # Final Velocity Objectives
        # self.d.Obj(self.final*1e3*(self.vz-vf[1])**2)
        # self.d.Obj(self.final*1e3*(self.vx-vf[0])**2)

        # Minimum Time Objective
        self.d.Minimize(1e4*self.tf)

        #solve
        # self.d.solve('http://127.0.0.1') # Solve with local apmonitor server
        self.d.open_folder()
        self.d.solve(disp=True)

        self.ts = np.multiply(self.d.time, self.tf.value[0])

        return self.a, self.u_p, self.ts

    def getTrajectoryData(self):
        return [self.ts, self.sx, self.sz, self.vx, self.vz, self.pitch, self.pitch_dot]

    def getInputData(self):
        return [self.ts, self.a]

# Main Code

opt = Optimizer()

s_ti = [0,0]
v_ti = 0
s_tf = [1000,500]
v_tf = [00.00, 00.0]
r_ti = 0 # inital orientation of the car
omega_ti = 0.0 # initial angular velocity of car

acceleration, turning, t_star = opt.optimize2D(s_ti, s_tf, v_ti, v_tf, r_ti, omega_ti)

# Printing stuff
# print('u', acceleration.value)
# print('tf', opt.tf.value)
# print('tf', opt.tf.value[0])
# print('u jump', opt.jump)
# for i in opt.u_jump: print(i.value)
print('u_jump', opt.u_jump.value)
print('jump his', opt.jump_hist.value)
print('v_mag', opt.v_mag.value)
print('a', opt.a.value)

# Plotting stuff

ts = opt.d.time * opt.tf.value[0]
t_max = opt.tf.value[0]
x_max = np.max(opt.sx.value)
vx_max = np.max(opt.vx.value)
z_max = np.max(opt.sz.value)
vz_max = np.max(opt.vz.value)
# plot results
fig = plt.figure(2)
ax = fig.add_subplot(111, projection='3d')
# plt.subplot(2, 1, 1)
Axes3D.plot(ax, opt.sx.value, ts, opt.sz.value, c='r', marker ='o')
plt.ylim(0, t_max)
plt.xlim(0, x_max)
plt.ylabel('time')
plt.xlabel('Position x')
ax.set_zlabel('position z')

n=5 #num plots
fig = plt.figure(3)
ax = fig.add_subplot(111, projection='3d')
# plt.subplot(2, 1, 1)
Axes3D.plot(ax, opt.vx.value, ts, opt.vz.value,  c='r', marker ='o')
plt.ylim(0, t_max)
plt.xlim(-1*vx_max, vx_max)
# plt.zlim(0, 2000)
plt.ylabel('time')
plt.xlabel('Velocity x')
ax.set_zlabel('vz')

plt.figure(1)
plt.subplot(n,1,1)
plt.plot(ts, opt.a, 'r-')
plt.ylabel('acceleration')

plt.subplot(n,1,2)
plt.plot(ts, np.multiply(opt.pitch, 1/math.pi), 'r-')
plt.ylabel('pitch orientation')

plt.subplot(n, 1, 3)
plt.plot(ts, opt.v_mag, 'b-')
plt.ylabel('vmag')

plt.subplot(n, 1, 4)
plt.plot(ts, opt.u_p, 'b-')
plt.ylabel('u_p')

plt.subplot(n, 1, 5)
plt.plot(ts, opt.u_jump, 'b-')
plt.plot(ts, opt.jump_hist, 'r-')
plt.ylabel('jump(b), jump hist(r)')

plt.show()

print('asdf')
 

Одна вещь, которую нужно попробовать, это решить с помощью IPOPT для инициализации, а затем APOPT, чтобы получить целочисленное решение. Еще одна вещь, которую нужно попробовать, это использовать MPCC для условия переключения, которое не зависит от двоичной переменной. Я обнаружил, что форма MPCC гораздо менее надежна, чем условие переключения двоичных переменных, потому что решатель часто застревает в седловой точке. Однако решение целых решений часто занимает гораздо больше времени.

 self.d.options.SOLVER=3
self.d.solve(disp=True)
self.d.options.TIME_SHIFT=0
self.d.options.SOLVER=1
self.d.solve(disp=True)
 

Вот решение с IPOPT:

 EXIT: Optimal Solution Found.

 The solution was found.

 The final value of the objective function is  506284.8987787149
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :  7.4613000000000005 sec
 Objective      :  506284.8987787149
 Successful solution
 ---------------------------------------------------
 

Целочисленное решение получается с помощью APOPT.

  --------- APM Model Size ------------
 Variable time shift OFF
 Number of state variables:    1286
 Number of total equations: -  1180
 Number of slack variables: -  40
 ---------------------------------------
 Degrees of freedom       :    66
 
 ----------------------------------------------
 Dynamic Control with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      2.72 NLPi:   92 Dpth:    0 Lvs:    3 Obj:  5.07E 05 Gap:       NaN
Iter:     2 I: -1 Tm:      0.53 NLPi:   17 Dpth:    1 Lvs:    2 Obj:  5.07E 05 Gap:       NaN
Iter:     3 I: -9 Tm:     47.59 NLPi:  801 Dpth:    1 Lvs:    1 Obj:  5.07E 05 Gap:       NaN
Iter:     4 I:  0 Tm:      2.26 NLPi:   35 Dpth:    1 Lvs:    3 Obj:  5.08E 05 Gap:       NaN
--Integer Solution:   2.54E 07 Lowest Leaf:   5.08E 05 Gap:   1.92E 00
Iter:     5 I:  0 Tm:      3.56 NLPi:   32 Dpth:    2 Lvs:    2 Obj:  2.54E 07 Gap:  1.92E 00
Iter:     6 I: -9 Tm:     54.65 NLPi:  801 Dpth:    2 Lvs:    1 Obj:  5.08E 05 Gap:  1.92E 00
Iter:     7 I: -1 Tm:      2.18 NLPi:   83 Dpth:    2 Lvs:    0 Obj:  5.08E 05 Gap:  1.92E 00
 No additional trial points, returning the best integer solution
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  113.5842 sec
 Objective      :  2.5419931399165962E 7
 Successful solution
 ---------------------------------------------------
 

APOPT предпочитает не переходить, чтобы минимизировать конечную цель. Возможно, вам потребуется добавить жесткое ограничение, vsum() равное of u_jump 1 . В курсе оптимизации есть дополнительные учебные пособия по MPCC и целочисленным / двоичным формам условий переключения.

Спасибо, что поделились своим приложением и держите нас в курсе!