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

**URL:** <https://discourse.julialang.org/t/ex-post-variable-values-not-equal-to-solver-computed-ones/2927>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [March 28, 2017, 2:13pm UTC](https://discourse.julialang.org/t/ex-post-variable-values-not-equal-to-solver-computed-ones/2927 "2017-03-28T14:13:51Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [March 28, 2017, 2:13pm UTC](https://discourse.julialang.org/t/ex-post-variable-values-not-equal-to-solver-computed-ones/2927/1 "2017-03-28T14:13:51Z")

</div>

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:

```julia
# 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:

```julia
------
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)`…