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

Question

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

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

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

# -*- 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()

标签： pythonoptimizationgekko

Answer 1

尝试使用m.Param()(or 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()