Fixed variable isn't fixed

[njgallagher]

I am trying to fix a model variable to run a base case scenario, but the value still varies. The variable I am fixing is M at parameter m0. Is there a different way I should be fixing the variable? Thank you.

The following gives:

julia> m0
Dict{String,Float64} with 2 entries:
"ercot" => 0.4
"ferc" => 0.8

julia> is_fixed.(M["ercot"])
true

julia> is_fixed.(M["ferc"])
true

I use the following code to fix:

for r in regions
fix(M[r], m0[r]; force = true)
end

After running I get:

julia> @show result_value.(M)
result_value.(M) = 1-dimensional DenseAxisArray{Float64,1,...} with index sets:
Dimension 1, ["ercot", "ferc"]
And data, a 2-element Array{Float64,1}:
133.48558627710796
76.66676884788406
1-dimensional DenseAxisArray{Float64,1,...} with index sets:
Dimension 1, ["ercot", "ferc"]
And data, a 2-element Array{Float64,1}:
133.48558627710796
76.66676884788406

julia> @show status
status = :StationaryPointFound
:StationaryPointFound

Here is the full code:

using Complementarity, JuMP

############# Parameters ####################
regions = ["ercot", "ferc"]
states = ["summer", "winter", "vortex"] # states of nature

σ = 0.25 # degree of relative risk aversion
esub = 0.25 # elasticity of subsitution between electricity and other goods
c0 = 31 # reference aggregate expenditures - electricity plus other goods

days = [300, 64, 1] # state occurance - days per year
πs = Dict(zip(states, days./365)) # state probabilities

maint = [0.4 0.8] # benchmark maintenance level
m0 = Dict(zip(regions, maint))

θ = 1/c0 # energy value share

scale_factor = [0 0] # maintenance cost scale factor
α = Dict(zip(regions, scale_factor))
γ = 0

###################### Shock Parameters ######################
d_shock = [ 1.0 1.2 1.2 ; # electricity demand shocks
1.0 1.0 1.5 ]

λ = Dict()
for i in 1:length(regions), j in 1:length(states) # demand shock table
λ[regions[i], states[j]] = d_shock[i,j]
end

w_shock = [ 0.0 0.2 0.25 ; # weather shocks
0.0 0.0 0.25 ]

w = Dict()
for i in 1:length(regions), j in 1:length(states) # weather shock table
w[regions[i], states[j]] = w_shock[i,j]
end
################### The Model ##########################

network = MCPModel() # declare model

################# Variables #################
@variable(network, K[r in regions] >= 0) # capacity
@variable(network, M[r in regions] >= 0) # maintenance

@variable(network, PE[r in regions, s in states] >= 0) # wholesale price
@variable(network, PR[r in regions] >= 0) # electricity generation resource

@variable(network, X[r in regions, s in states] >= 0) # sales to other regions
@variable(network, PX[s in states]) # traded electricity price

#################### Expressions #############################
@NLexpression(network, PC[r in regions, s in states], # state-contingent price
(θ * (PE[r,s] * λ[r,s])^(1-esub) + 1-θ)^(1/(1-esub)))
@NLexpression(network, PEU[r in regions], #
sum(πs[s] * PC[r,s]^(1-σ) for s in states)^(1/(1-σ)))
@NLexpression(network, EU[r in regions],
1/PEU[r])
@NLexpression(network, C[r in regions, s in states],
EU[r] * (PEU[r]/PC[r,s])^σ)
@NLexpression(network, CM[r in regions],
K[r] * α[r]/(1-γ) * (1-(1-M[r])^(1-γ)))

#################### Equations #############################
@mapping(network, profit_K[r in regions],
sqrt(PR[r]) - sum(πs[s] * PE[r,s] * (1-w[r,s](1-M[r])) for s in states))
@mapping(network, profit_M[r in regions],
α[r]/(1-M[r])^γ - sum(πs[s] * PE[r,s] * w[r,s] for s in states))
@mapping(network, profit_X[r in regions, s in states],
PE[r,s] - PX[s])
@mapping(network, market_PE[r in regions, s in states],
K[r] * (1-w[r,s](1-M[r])) -
C[r,s] * λ[r,s] * (PC[r,s]/(PE[r,s]*λ[r,s]))^esub + X[r,s])
@mapping(network, market_PR[r in regions],
0.5 * (1 - K[r]/sqrt(PR[r])))
@mapping(network, market_PX[s in states],
sum(X[r,s] for r in regions))

############## Assign Complementarity #################
@complementarity(network, profit_K, K)
@complementarity(network, profit_M, M)
@complementarity(network, profit_X, X)
@complementarity(network, market_PE, PE)
@complementarity(network, market_PR, PR)
@complementarity(network, market_PX, PX)

############## Starting Variable Values######################
for r in regions
set_start_value(K[r], 1.0)
set_start_value(PR[r], 1.0)
fix(M[r], m0[r]; force = true)
for s in states
set_start_value(PE[r,s], 1.0)
end
end

status = solveMCP(network; convergence_tolerance=1e-8, output="yes", time_limit=3600)

[jonmbecker]

I noticed this as well, not sure what I'm doing wrong. I ended up setting upper and lower bounds in variable declaration as an interim solution. Alternatively, you should be able to set the upper and lower bound separately after declaration.

@variable(model,1.0>=PFX>=1.0, start=1)