Почему этот сценарий GEKKO не дает лучшего решения?

Вот мой код. Я максимизирую выражение abs (expr1) с учетом ограничения abs (expr1) = abs (expr2).

import numpy as np
from gekko import GEKKO

#init
m = GEKKO(remote=False)

x2,x3,x4,x5,x6,x7,x8 = [m.Var(lb=-2*np.pi, ub=2*np.pi) for i in range(7)]

#constraint
m.Equation((1/8)*m.abs(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*m.sin((1/2)*x7)*((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*m.sin(x3)*m.sin(x8))+m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*m.cos(x5)*m.cos(x8)+(-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin(x3)*m.sin(x8)))) == (1/8)*m.abs(m.sin((1/2)*x6)*(4*m.cos(x3)*m.cos(x8)*(m.cos((1/2)*x5)*m.cos((1/2)*x7)*m.sin(x2)*m.sin(x4)+2*m.cos(x2)*m.cos(x5)*m.sin((1/2)*x4)*m.sin((1/2)*x7))+(-1)*m.sin(x3)*((3+m.cos(x5))*m.cos((1/2)*x7)*m.sin(x2)*m.sin(x4)+m.cos(x2)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin((1/2)*x4)*m.sin((1/2)*x7))*m.sin(x8))))

#objective
m.Obj(-((1/8)*m.abs(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*m.sin((1/2)*x7)*((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*m.sin(x3)*m.sin(x8))+m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*m.cos(x5)*m.cos(x8)+(-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin(x3)*m.sin(x8))))))

#Set global options
m.options.IMODE = 3 

#execute
m.solve()

#output
print('')
print('Results')
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
print('x5: ' + str(x5.value))
print('x6: ' + str(x6.value))
print('x7: ' + str(x7.value))
print('x8: ' + str(x8.value))

Решение, которое я получаю, - это все x_i = 0, что является допустимым, но не лучшим решением. Например

x2,x3,x4,x5,x6,x7,x8 = 0.9046, 1.9540, 1.8090, 0, 1.8090, 6.2832, 4.3291

удовлетворяет ограничению (до 5-го знака после запятой), а цели достигают -0,3003 (что все еще не лучший результат, но это пример). Я пытался поиграть с вариантами допуска, но безуспешно. Обратите внимание, что если я удалю ограничение равенства, цель должным образом достигнет максимального значения -1.

Почему решатель застревает на нулевом решении?


gekko nonlinear-optimization
person optiman    schedule 06.10.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Решатели в Gekko - это локальные минимизаторы, а не глобальные минимизаторы. У вашей задачи много локальных минимумов с функциями sin и cos. Вы можете получить глобальный минимум, используя метод нескольких запусков или глобальный подход, такой как смоделированный отжиг. Я использовал модификацию вашего скрипта с многозапуском. Случайные значения от -2*np.pi до 2*np.pi используются для инициализации переменных.

for xi in x:
    xi.value = np.random.rand(20)*4*np.pi - 2*np.pi

IMODE=2 решает все эти случаи одновременно.

m.options.IMODE = 2

Если вам нужно выполнить много случаев, вы можете распараллелить этот расчет с несколькими потоками . Вам также следует переключиться на m.abs3 вместо m.abs, чтобы избежать проблем с прерывистой производной на нуле. Другая стратегия - возвести обе части уравнения в квадрат, чтобы избежать абсолютного значения. Вот полная версия:

import numpy as np
from gekko import GEKKO

#init
m = GEKKO(remote=False)

x = [m.Var(lb=-2*np.pi, ub=2*np.pi) for i in range(7)]
for xi in x:
    xi.value = np.random.rand(20)*4*np.pi - 2*np.pi
x2,x3,x4,x5,x6,x7,x8 = x

#constraint
m.Equation((1/8)*m.abs3(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*m.sin((1/2)*x7)*\
    ((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*m.sin(x3)*m.sin(x8))+\
    m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*m.cos(x5)*m.cos(x8)+\
    (-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin(x3)*m.sin(x8)))) == \
    (1/8)*m.abs3(m.sin((1/2)*x6)*(4*m.cos(x3)*m.cos(x8)*(m.cos((1/2)*x5)*\
    m.cos((1/2)*x7)*m.sin(x2)*m.sin(x4)+2*m.cos(x2)*m.cos(x5)*m.sin((1/2)*\
    x4)*m.sin((1/2)*x7))+(-1)*m.sin(x3)*((3+m.cos(x5))*m.cos((1/2)*x7)*\
    m.sin(x2)*m.sin(x4)+m.cos(x2)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*\
    m.sin((1/2)*x4)*m.sin((1/2)*x7))*m.sin(x8))))

#objective
obj = m.Intermediate(-((1/8)*m.abs3(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*\
    m.sin((1/2)*x7)*((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*\
    m.sin(x3)*m.sin(x8))+m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*\
    m.cos(x5)*m.cos(x8)+(-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*\
    m.sin(x3)*m.sin(x8))))))
m.Obj(obj)

#Set global options
m.options.IMODE = 2

#execute
m.solve()

#output
print('')
print('Best Result')
i = np.argmin(obj.value)
print('x2: ' + str(x2.value[i]))
print('x3: ' + str(x3.value[i]))
print('x4: ' + str(x4.value[i]))
print('x5: ' + str(x5.value[i]))
print('x6: ' + str(x6.value[i]))
print('x7: ' + str(x7.value[i]))
print('x8: ' + str(x8.value[i]))
print('obj: ' + str(obj.value[i]))

Есть несколько лучших решений с целью -0,5.

Best Result
x2: -3.1415936876
x3: 6.2545093655
x4: 3.1415896007
x5: -2.0848973806e-05
x6: -4.7122128433
x7: -4.712565114
x8: 0.029076147797
obj: -0.50000008426

Best Result
x2: -3.1416640191
x3: 3.1415941185
x4: 3.1415948958
x5: -3.1416088732
x6: 1.5708701192
x7: -4.7124627728
x8: -3.1415893349
obj: -0.5000000992
person John Hedengren    schedule 08.10.2020
comment
Влияет ли квадрат с обеих сторон на скорость? Что делать, если ограничение сложное и мне нужно «сопрягать» квадрат? Он все еще быстрее абс3? - person optiman; 08.10.2020
comment
Квадраты членов делают модель более нелинейной. Функция abs3 добавляет двоичную переменную, которая может быть сложной для крупномасштабных систем уравнений. Я рекомендую вам попробовать оба подхода. - person John Hedengren; 09.10.2020