Ex-post variable values not equal to solver-computed ones

jump

#1

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)