Python GEKKO：如何在 gekko 中表达随机变量的约束？

Question

我花了一些时间尝试用 gekko 解决以下问题，但很难判断我是否正确编写了代码。我这么说是因为当我稍微减少 tau 的值时，问题就变得不可行（理论上并非如此）。我需要你的帮助；

这是一个跟踪问题，状态变量取决于多个值 xi。

这是我的代码；

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

# Initialize model
m = GEKKO(remote=False)

# Time points
n = 10
m.time = np.linspace(0, 10, n)

# parameters
alpha = 0.95
m1 = 0.5
m2 = 0.3
tau = 0.5

# Number of random samples
N = 10

# create the samples
Xi = np.random.rand(N,n)
S = []
for i in range(N):
  S.append(m.Param(list(Xi[i])))

# State variables
x = m.Array(m.Var, N, value=5.0)

# Control variable
u = m.Var()

# reference trajectory
x_val = [5.0, 2.31, 1.45, 1.01, 0.75, 0.57, 0.45, 0.36, 0.31, 0.28]
x_r = m.Param(x_val)

# Equation
m.Equations(5 * x[i].dt() == -x[i]**2 + u + S[i] for i in range(N)) 

m.Equation((1/N)*m.sum([(1 + m1 * tau) / (1 + m2 * tau * m.exp((-(x[i] - x_r[i])**2 + 10) / tau)) for i in range(N)]) <= 1 - alpha)

m.Obj((1/N)*(m.sum([(x[i] - x_r[i])**2 for i in range(N)])) + u**2)

# Solve the differential equation
m.options.IMODE = 6  #
m.options.SOLVER = 3; # 1: APOPT, 2:BPOPT, 3:IPOPT
m.options.MAX_ITER = 10000 # change maximum iterations
m.solve(disp=True)

我预计随着 tau 值的降低，目标值也会降低。

谢谢！

Answer 1

看起来图像中的方程和模型中写的方程有差异。

m.Equation((1/N)*m.sum([(1 + m1 * tau) / 
           (1 + m2 * tau * m.exp((-(x[i] - x_r[i])**2 + 10) /
           tau)) for i in range(N)]) <= 1 - alpha)

分母中多了一个

/tau

，符号应为

(xi-xd)^2-10

。

收敛的另一个潜在问题是被零除。提高收敛性的一种方法是将分母限制为大于一个小值，例如

1e-3

。

den = m.Array(m.Var,N,lb=1e-3)
m.Equations([den[i]==1+m2*tau*m.exp(((x[i]-x_r[i])**2-10)) for i in range(N)])
m.Equation((1/N)*m.sum([(1 + m1 * tau) / den[i] for i in range(N)]) <= 1 - alpha)

不可行可能是由于不平等约束。尝试用

p-alpha

和更大的

p

值放松不等式约束，直到问题变得可行。在这种情况下，当

p=5

时，就找到了成功的解决方案。

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

# Initialize model
m = GEKKO(remote=True)

# Time points
n = 10
m.time = np.linspace(0, 10, n)

# parameters
alpha = 0.95
m1 = 0.5
m2 = 0.3
tau = 0.5
p = 5.0

# Number of random samples
N = 10

# create the samples
Xi = np.random.rand(N,n)
S = []
for i in range(N):
  S.append(m.Param(list(Xi[i])))

# State variables
x = m.Array(m.Var, N, value=5.0)

# Control variable
u = m.Var()

# reference trajectory
x_val = [5.0, 2.31, 1.45, 1.01, 0.75, 0.57, 0.45, 0.36, 0.31, 0.28]
x_r = m.Param(x_val)

# Equation
m.Equations(5 * x[i].dt() == -x[i]**2 + u + S[i] for i in range(N))

den = m.Array(m.Var,N,lb=1e-3)
m.Equations([den[i]==1+m2*tau*m.exp(((x[i]-x_r[i])**2-10)) for i in range(N)])
m.Equation((1/N)*m.sum([(1 + m1 * tau) / den[i] for i in range(N)]) <= p - alpha)

m.Obj((1/N)*(m.sum([(x[i] - x_r[i])**2 for i in range(N)])) + u**2)

# Solve the differential equation
m.options.IMODE = 6  #
m.options.SOLVER = 3 # 1: APOPT, 2:BPOPT, 3:IPOPT
m.options.MAX_ITER = 200 # change maximum iterations
m.solve(disp=True)