Проблема объезда препятствий с помощью GEKKO (Решение не найдено, и решение существует)

Исходный вопрос можно найти здесь.

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

Обычный импорт:

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

Инициализация модели и определение моментов времени:

m = GEKKO() # initialize gekko
nt = 501
m.time = np.linspace(0.0,1.0,nt)
t = m.Param(value=m.time)

p = np.zeros(nt) # mark final time point
p[-1] = 1.0
final = m.Param(value=p)

Определение переменных состояния и управления:

x = m.Var(value=0.0)
y = m.Var(value=0.0)
z = m.Var(value=0.0)

V = m.FV(value=2.0, lb=0.0, ub=100.0)
V.STATUS = 1

theta = m.CV(1.0)
theta.STATUS = 1
theta.SPHI = 0
theta.SPLO = 6.28

V был выбран в качестве фиксированной переменной, поскольку в формулировке задачи упоминается, что:

V - постоянная скалярная скорость.

Предполагается, что theta находится в [0,2*pi], поскольку он не упоминается напрямую в формулировке проблемы.

z используется в качестве прокси-переменной для реализации интегральной цели.

Уравнения связи определены:

m.Equation(x.dt()==V*m.cos(theta))
m.Equation(y.dt()==V*m.sin(theta))
m.Equation(x*final==1.2)
m.Equation(y*final==1.6)
m.Equation(((x-0.4)**2)+((y-0.5)**2)>=0.1)
m.Equation(((x-0.8)**2)+((y-1.5)**2)>=0.1)

Определение переменных, необходимых для получения производных 2-го порядка:

dx = m.Var(value = 0.0)
dy = m.Var(value = 0.0)
ddx = m.Var(value = 0.0)
ddy = m.Var(value = 0.0)

Соответствующие уравнения для работы:

m.Equation(dx==x.dt())
m.Equation(dy==y.dt())
m.Equation(ddx==dx.dt())
m.Equation(ddy==dy.dt())

Установка производных прокси-переменной в качестве цели внутри интеграла и решение:

m.Equation(z.dt()==(ddx**2)+(ddy**2))
m.Obj(z*final)

m.options.IMODE = 6 # optimal control mode
m.solve(disp=True)

Я не совсем уверен, верен ли способ, которым я закодировал данный вопрос. Как видно по данной ссылке, решение проблемы существует. Я пробовал использовать APOPT и IPOPT для решения этой проблемы, но оба решателя не смогли найти решение. COLD_START нельзя использовать с DOF < 0. Любая помощь с этим будет принята с благодарностью!


gekko
person fakeshade    schedule 21.03.2021    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Проблема с моделью - это два конечных ограничения, потому что каждая точка, кроме конечной точки (где final = 1) - это 0=1.2 и 0=1.6.

m.Equation(x*final==1.2)
m.Equation(y*final==1.6)

Требуется альтернативная форма.

m.Equation((x-1.2)*final==0)
m.Equation((y-1.6)*final==0)

Ниже представлено возможное, но не оптимизированное решение.

планирование пути

import numpy as np
from gekko import GEKKO

m = GEKKO()
nt = 51
m.time = np.linspace(0,1,nt)

# mark final time point
p = np.zeros(nt); p[-1] = 1 
final = m.Param(p)

# define theta as a Variable
#theta = m.Var(value=np.pi/2.0,lb=0,ub=2*np.pi)

# define theta as an MV
theta = m.MV(value=0,lb=0,ub=2*np.pi)
theta.STATUS = 1; theta.DCOST = 0

V = m.FV(value=2.15,lb=0,ub=10) # velocity
#V.STATUS = 1

x,y,dx,dy,ddx,ddy = m.Array(m.Var,6)

m.Equation(x.dt()==V*m.cos(theta))
m.Equation(y.dt()==V*m.sin(theta))
m.Equation(((x-0.4)**2)+((y-0.5)**2)>=0.1)
m.Equation(((x-0.8)**2)+((y-1.5)**2)>=0.1)

m.Equation(dx==x.dt())
m.Equation(dy==y.dt())
m.Equation(ddx==dx.dt())
m.Equation(ddy==dy.dt())

# alternative 1 terminal constraint
m.Minimize(1e4*(x-1.2)**2)
m.Minimize(1e4*(y-1.6)**2)

# alternative 2 terminal constraint
#m.fix_final(x,1.2)
#m.fix_final(y,1.6)

# alternative 3 terminal constraint
#m.Equation((x-1.2)*final==0)
#m.Equation((y-1.6)*final==0)

# objective
m.Minimize(m.integral((ddx**2)+(ddy**2))*final)

m.options.SOLVER = 3
m.options.MAX_ITER = 1000
m.options.NODES = 3
m.options.IMODE = 6 # optimal control mode
m.solve(disp=True)

import matplotlib.pyplot as plt
plt.figure(figsize=(5,5))
th = np.linspace(0,2*np.pi,500)
x1 = np.sqrt(0.1)*np.cos(th)+0.4
y1 = np.sqrt(0.1)*np.sin(th)+0.5
x2 = np.sqrt(0.1)*np.cos(th)+0.8
y2 = np.sqrt(0.1)*np.sin(th)+1.5
plt.plot(x1,y1,'r',x2,y2,'r');
plt.plot(x.value,y.value,'b-o',markersize=2)
plt.plot([1.2],[1.6],'o',color='orange')
plt.xlabel('x'); plt.xlim([0,2.0])
plt.ylabel('y'); plt.ylim([0,2.0])
plt.title('Obstacle avoidance state variables, Speed = %2.4g'%V.value[0])

IPOPT не может решить проблему, в которой V.STATUS=1 (вычислено). Я попробовал еще пару ограничений терминала, но решатели (APOPT и IPOPT) не увенчались успехом.

person John Hedengren    schedule 26.03.2021