Complementarity.jl icon indicating copy to clipboard operation
Complementarity.jl copied to clipboard

Fixed variable isn't fixed

Open njgallagher opened this issue 3 years ago • 1 comments

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)

njgallagher avatar Mar 19 '21 17:03 njgallagher

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)

jonmbecker avatar Jun 08 '21 00:06 jonmbecker