I have a simple JuMP optimisation problem where the objective is to maximise forest harvesting under a budget constraint of subsides/tax to push up/down harvesting, a supply curve that depends on price+subsides and forest stock levels and a demand curve that depends from price.
The problem is that when I compute demand and supply ex-post (from price, forest stocks, subsides and parameters) I get very different results than the demand and supply reported by the solver/jump.
Unless I done something idiot in the code (I checked several times, but…) it seems that JuMP is ignoring this part of the problem…
This is the code:
# Import of JuMP
using JuMP, Ipopt
#####################################
# Sets
#####################################
time = collect(1:100)
#####################################
# Parameters
#####################################
# Forest
# ------
# Parameters of the logistic function S = m / (1+n*e^-lt) -> G = lS - (l/m) S^2
l = 0.05 # growth rate
m = 8000. # maximum value
n = 600. # ?
o = l/m
SZ =2500. # initial volumes M m^3
# Markets
# -------
# Parameters of the supply function R = a * (p*x)^α * S^σ
a = 0.8 # coefficient
alpha = 0.2 # price exponent
sigma = 0.5 # volumes exponent
# Parameters of the demand function D = b * p^β
b = 400. # coefficient
beta = -0.5 # price exponent
Bg = SZ * 30. * length(time) # budget M €
r = 0.00 # mkt interest rate
rho = 0.00 # social discount rate
#####################################
# Model declaration
#####################################
m = Model(solver=IpoptSolver(print_level=0)) #print_level=0
#####################################
# Variables
#####################################
@variable(m, 0.1 <= G[t in time[2:end]] <= 200, start=110. ) # growth at time t
@variable(m, 1.0 <= S[t in time] <= 15000, start=2500. ) # stocks at time t
@variable(m, 0.1 <= R[t in time[2:end]] <= 500, start=77. ) # harvesting at time t
@variable(m, x[t in time[2:end]], start=1.0 ) # policy subside rate
@variable(m, 0.1 <= p[t in time[2:end]] <= 300, start=30. ) # price on each year
@variable(m, 0.1 <= D[t in time[2:end]] <= 500, start=77. ) # demand each year
#####################################
# Constraints
#####################################
@NLconstraint(m, growth[t in time[2:end]], l * S[t-1] - o * (S[t-1])^2 == G[t] ) # natural growth
@NLconstraint(m, harvesting[t in time[2:end]], a * (p[t] * x[t])^alpha * (S[t-1]^sigma) == R[t]) # harvesting level ≡ supply
@NLconstraint(m, demand[t in time[2:end]], b * p[t]^beta == D[t]) # demand
@NLconstraint(m, mktclearance[t in time[2:end]], D[t] == R[t]) # mkt (static) clearance
@NLconstraint(m, budget, sum((x[t]-1)*p[t]/(1+r)^t for t in time[2:end]) <= Bg) # financial public budget
@NLconstraint(m, stocksVar[t in time[2:end]], S[t] == S[t-1] + G[t] - R[t]) # overall stock variation (stocks are at the end of the year)
@NLconstraint(m, finalStocks, S[time[1]] <= S[time[end]]) # final stocks
@NLconstraint(m, initialStocks, S[time[1]] == SZ) # initial stocks
#####################################
# Objective
#####################################
@NLobjective(m, :Max, sum(R[t]/(1+rho)^t for t in time[2:end]))
#print(m)
#####################################
# Solve
#####################################
status = solve(m)
#####################################
# Results
#####################################
if status == :Optimal
x_out = getValue(x)
p_out = getValue(x)
R_out = getValue(R)
D_out = getValue(D)
S_out = getValue(S)
println("------")
println("p: "*string(p_out[50]))
println("x: "*string(x_out[50]))
println("S: "*string(S_out[49]))
println("a: "*string(a))
println("alpha: "*string(alpha))
println("sigma: "*string(sigma))
println("b: "*string(b))
println("beta: "*string(beta))
println("R: "*string(R_out[50]))
println("D: "*string(D_out[50]))
R_comp = a * (S_out[49]^sigma) * ((p_out[50] * x_out[50])^alpha)
D_comp = b * (p_out[50]^beta)
println("R_comp: "*string(R_comp))
println("D_comp: "*string(D_comp))
else
println("Model didn't solved")
println(status)
end
and this is the output:
------
p: 1.884864366117979
x: 1.884864366117979
S: 4000.000000031871
a: 0.8
alpha: 0.2
sigma: 0.5
b: 400.0
beta: -0.5
R: 100.00000000000476
D: 100.00000000000476
R_comp: 65.19766232244544
D_comp: 291.3532985436304
As you can see, both reported R (harvesting == supply) and D are different than computed ones… Any guess on the reason ??
EDITED: it was actually a typo: p_out = getValue(p) instead of p_out = getValue(x)…