如何在Gekko中编写累积密度函数和期望值？

Question

假设

X

是

[\mu-a,\mu+a]

范围内的均匀随机变量。如果

a=0

，那么

X

将是一个常数。此外，

X

的支持可能是正数，也可能是负数，或者既是负数又是正数。需要注意的是，

\mu

可以是负值也可以是正值。假设在我们的优化问题中，决策变量是

y>= 0

。现在，我的目标函数如下：

P = (c_1 + F_X(y))*(c_2+E(X I{X>=y}))

其中

c_i

是常数，

E(X)

是期望值，

I{}

是指标函数。我想知道如何优化这个问题。另外，

F_X(y)

是随机变量

X

的CDF。

Answer 1

有多种方法可以解决这个问题，具体取决于是否有数据点或者这是否是一个严格的理论问题。

理论优化

对于数学问题，很容易将目标函数可视化为单个可调变量

y

的函数。如果

y<=mu-a

则 CDF 函数为零，如果

y>=mu+a

则 CDF 函数为 1。预期值

E(X 1{X>=y})

包括指示函数

1{X>=y}

，如果

X>=y

则为 1，否则为零。该函数在 0 到

E(X 1{X>=y}) = 1/(4a) * ((u+a)^2-y^2)

范围内变为

mu

。该图显示最大目标约为

y=4.2

。

import numpy as np
import matplotlib.pyplot as plt

# Parameters
mu = 5; a = 2
c_1 = 1.0; c_2 = 2.0

y = np.linspace(0,10,200)
F_X = np.clip(1/(2*a) * (y-(mu-a)),0,1)
E = np.clip((1/(4*a))*((mu+a)**2-y**2),0,mu)

P = (c_1+F_X)*(c_2+E)

plt.figure(figsize=(6,4))
plt.subplot(2,1,1)
plt.plot(y,P,'k-.',label='P')
plt.grid(); plt.legend()
plt.ylabel('Objective')
plt.subplot(2,1,2)
plt.plot(y,F_X,'b:',label='F_X (CDF)')
plt.plot(y,E,'r--',label='E (Expected Value X>=y)')
plt.xlabel('y'); plt.legend(); plt.grid()
plt.tight_layout()
plt.savefig('obj.png',dpi=300)
plt.show()

Gekko 可用于精确求函数的最大值。

from gekko import GEKKO

# Parameters
mu = 5; a = 2
c_1 = 1.0; c_2 = 2.0

m = GEKKO()

# variable with bound, y>=0
y = m.Var(mu,lb=0)

# CDF
F_X = m.Var(lb=0,ub=1)
m.Equation(F_X == (1/(2*a))*(y-(mu-a)))

# E(X 1{X>=y})
E = m.Var(lb=0,ub=mu)
m.Equation(E == (1/(4*a))*((mu+a)**2-y**2))

m.Maximize((c_1+F_X)*(c_2+E))
m.solve(disp=False)

print(f'y={y.value[0]}')
print(f'F_X={F_X.value[0]}')
print(f'E={E.value[0]}')

解决办法是：

y=4.3333337334
F_X=0.33333343336
E=3.7777774075

数据驱动优化

对于数据点，使用代理函数来近似 CDF。如果

X

均匀分布，它会接近一条具有足够数据的线：

F_X(y) = 1/(2a) * (y - (mu-a))

但是，如果数据分布不均匀或者没有足够的数据来确定

F_X(y)

，则可以凭经验确定 CDF 并使用分段线性 (PWL) 函数进行近似，例如：

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

# Parameters
mu = 5
a = 2
num_points = 100

# Generate 20 data points for X
# uniform random variable within [mu-a, mu+a]
a = abs(a) # if "a" is negative
np.random.seed(14)
if a == 0:
    X = np.full(num_points, mu)
else:
    X = np.random.uniform(mu-a, mu+a, num_points)

# Sort data points
X_sorted = np.sort(X)

# Generate CDF values
Z_cdf = np.arange(1, num_points + 1) / num_points

# Model
m = GEKKO()
m.options.IMODE = 2

# Piece-wise linear interpolation
x_pwl = m.Param(value=np.linspace(min(X_sorted), max(X_sorted), 100))
z_pwl = m.Var(lb=0,ub=1)
m.pwl(x_pwl, z_pwl, X_sorted, Z_cdf)

m.solve(disp=False)

plt.figure(figsize=(6, 3.5))
plt.plot(X_sorted, Z_cdf, 'bo', label='Data Points')
plt.plot([mu-a,mu+a], [0,1], 'k:', label='CDF Limit with Large Uniform Data')
plt.plot(x_pwl.value, z_pwl.value, 'r--', label='Piecewise Linear')
plt.xlabel('X'); plt.ylabel('CDF'); plt.legend(); plt.grid()
plt.tight_layout(); plt.savefig('pwl.png',dpi=300)
plt.show()

一旦有了 CDF 的函数逼近（线性或 PWL），就可以用来优化函数。

m.if3()

gekko 函数处理条件语句。但是，不需要此条件函数，因为最大值始终在下限

mu-a

和上限

mu+a

的限制内。