python - 在动态优化问题中实现不确定性

Question

我有一个关于在 Gekko 优化问题中实现不确定性项的问题。我是编码初学者，从添加“鱼类管理”示例中的部分开始。我有两个主要问题。

1. 我想在模型中添加一个不确定性项(例如， future 价格波动)，但我似乎不了解该模块的工作原理。我正在尝试从某个分布中抽取一个随机值并将其放入 m.Var，'ss'，希望随着 t 的移动，模块每次都会获取每个值。但似乎该模块不能那样工作。我想知道是否有任何方法可以在流程中实现不确定性条款。

2. 假设在使用 A 和 E 之间分配初始土地 A(0) 的优化问题通过控制土地转化为单个代理解决，e，我计划将其扩展到多个代理问题。例如，如果代理之间的土地质量 h 和土地数量 A 不同，n，我计划通过调用初始 m.Var 值和来自 a 的一些参数来解决使用算法的多个优化问题加载的数据框。如果可以的话，我可以对这个计划做一个简短的评论吗？

# -*- coding: utf-8 -*-

from gekko import GEKKO
from scipy.stats import norm
from scipy.stats import truncnorm
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import operator
import math
import random

# create GEKKO model
m = GEKKO()


# Below, an agent is given initial land A(0) and makes a decision to convert this land to E
# Objective of an agent is to get maximum present utility(=log(income)) from both land use, income each period=(A+E(1-y))*Pa-C*u+E*Pe
# Uncertainty in future price lies for Pe, which I want to include with ss below
# After solving for a single agent, I want to solve this for all agents with different land quality h, risk aversion, mu, and land size A
# Then lastly collect data for total land use over time


# time points
n=51
year=10
k=50
m.time = np.linspace(0,year,n)
t=m.time
tt=t*(n-1)/year
tt = tt.astype(int)
ttt = np.exp(-t/(n-1))

# constants
Pa = 1
Pe = 1
C = 0.1
r = 0.05
y = 0.1

# distribution
# I am trying to generate a distribution, and use it as uncertainty term later in objective function
mu, sigma = 1, 0.1 # mean and standard deviation
#mu, sigman = df.loc[tt][2]
sn = np.random.normal(mu, sigma, n)
s= pd.DataFrame(sn)
ss=s.loc[tt][0]

# Control
# What is the difference between MV and CV? They give completely different solution
# MV seems to give correct answer
u = m.MV(value=0,lb=0,ub=10)
u.STATUS = 1
u.DCOST = 0
#u = m.CV(value=0,lb=0,ub=10)

# Variables
# m.Var and m.SV does not seem to lead to different results
# Can I call initial value from a dataset? for example, df.loc[tt][0] instead of 10 below? 
# A = m.Var(value=df.loc[tt][0])
# h = m.Var(value=df.loc[tt][1])

A = m.SV(value=10)
E = m.SV(value=0)
#A = m.Var(value=10)
#E = m.Var(value=0)
t = m.Param(value=m.time)
Pe = m.Var(value=Pe)
d = m.Var(value=1)

# Equation
# It seems necessary to include restriction on u
m.Equation(A.dt()==-u)
m.Equation(E.dt()==u)
m.Equation(Pe.dt()==-1/k*Pe)
m.Equation(d==m.exp(-t*r))
m.Equation(A>=0)


# Objective (Utility)
J = m.Var(value=0)

# Final objective
# I want to include ss, uncertainty term in objective function
Jf = m.FV()
Jf.STATUS = 1
m.Connection(Jf,J,pos2='end')
#m.Equation(J.dt() == m.log(A*Pa-C*u+E*Pe))
m.Equation(J.dt() == m.log((A+E*(1-y))*Pa-C*u+E*Pe)*d)
#m.Equation(J.dt() == m.log(A*Pa-C*u+E*Pe*ss)*d)

# maximize profit
m.Maximize(Jf)
#m.Obj(-Jf)

# options
m.options.IMODE = 6  # optimal control
m.options.NODES = 3  # collocation nodes
m.options.SOLVER = 3 # solver (IPOPT)

# solve optimization problem
m.solve()

# print profit
print('Optimal Profit: ' + str(Jf.value[0]))


# collect data
# et=u.value
# print(et)
# At=A.value
# print(At)
# index = range(1, 2)
# columns = range(1, n+1)
# Ato=pd.DataFrame(index=index,columns=columns)
# Ato=At

# plot results
plt.figure(1)
plt.subplot(4,1,1)
plt.plot(m.time,J.value,'r--',label='profit')
plt.plot(m.time[-1],Jf.value[0],'ro',markersize=10,\
         label='final profit = '+str(Jf.value[0]))
plt.plot(m.time,A.value,'b-',label='Agricultural Land')
plt.ylabel('Value')
plt.legend()
plt.subplot(4,1,2)
plt.plot(m.time,u.value,'k.-',label='adoption')
plt.ylabel('conversion')
plt.xlabel('Time (yr)')
plt.legend()
plt.subplot(4,1,3)
plt.plot(m.time,Pe.value,'k.-',label='Pe')
plt.ylabel('price')
plt.xlabel('Time (yr)')
plt.legend()
plt.subplot(4,1,4)
plt.plot(m.time,d.value,'k.-',label='d')
plt.ylabel('Discount factor')
plt.xlabel('Time (yr)')
plt.legend()
plt.show()

Answer 1

尝试使用m.Param()(或m.FV、m.MV)代替m.Var () 实现为每个优化指定的不确定值。使用 m.Var()，可以指定初始值，但随后由优化器计算初始值，初始猜测值不再保留。

另一种实现不确定性的方法是在实例之间创建一个具有不同参数的模型数组。这是一个简单的示例，其中 K 作为随机选择的值。这导致模型的 10 个实例被优化为值 40。

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

# uncertain parameter
n = 10
K = np.random.rand(n)+1.0

m = GEKKO()
m.time = np.linspace(0,20,41)

# manipulated variable
p = m.MV(value=0, lb=0, ub=100)
p.STATUS = 1 
p.DCOST = 0.1  
p.DMAX = 20

# controlled variable
v = m.Array(m.CV,n)
for i in range(n):
    v[i].STATUS = 1
    v[i].SP = 40
    v[i].TAU = 5
    m.Equation(10*v[i].dt() == -v[i] + K[i]*p)

# solve optimal control problem
m.options.IMODE = 6
m.options.CV_TYPE = 2
m.solve()

# plot results
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',LineWidth=2)
plt.ylabel('MV')
plt.subplot(2,1,2)
plt.plot([0,m.time[-1]],[40,40],'k-',LineWidth=3)
for i in range(n):
    plt.plot(m.time,v[i].value,':',LineWidth=2)
plt.ylabel('CV')
plt.xlabel('Time')
plt.show()