Implementing uncertainty in dynamic optimization problem

user_15312117 0
Viewed 84

I have a question about implementing uncertainty terms into the Gekko optimization problem. I am a beginner in coding and starting by adding parts from the "fish management" example. I have two main questions.

  1. I would like to add an uncertainty term in the model (e.g. fluctuating future prices) but it seems like I am not understanding how the module works. I am trying to draw a random value from a certain distribution and put it into m.Var, 'ss', hoping the module will take each value at each time as t moves on. But it seems like the module does not work that way. I am wondering if there is any way I can implement uncertainty terms into the process.

  2. Assuming the optimization problem to allocate initial land, A(0), between use A and E, is solved for a single agent by controlling land to convert, e, I am planning to expand this to a multiple agents problem. For example, if land quality, h, and land quantity A differ among agents, n, I am planning to solve multiple optimization problems using for algorithm by calling the initial m.Var value and some parameters from a loaded dataframe. If possible, may I have a brief comment on this plan?

# -*- 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()
gekko optimization python
1 Answers

Try using the m.Param() (or m.FV, m.MV) instead of m.Var() to implement an uncertain value that is specified for each optimization. With m.Var(), an initial value can be specified but then it is calculated by the optimizer and the initial guess is no longer preserved.

Another way to implement uncertainty is to create an array of models with different parameters among the instances. Here is a simple example with K as a randomly selected value. This leads to 10 instances of the model that is optimized to a value of 40.

Uncertainty

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()
user_23669410
Related