chkwon / Complementarity.jl

provides a modeling interface for mixed complementarity problems (MCP) and math programs with equilibrium problems (MPEC) via JuMP
Other
75 stars 20 forks source link

Setting start values for Array Variables #58

Closed njgallagher closed 3 years ago

njgallagher commented 3 years ago

I want to set starting values for Array variables. Can I do this in one line, or does it need to happen individually? Below is the code and the Error report. I noted where the first error is. If I should post this issue differently please let me know.

I am also wondering if I can set starting value for expressions. See lines 101-105.

############ The File #####################
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 ]

d_shock1 = [ 1.0 1.0 1.0 ;   # no electricity demand shocks - base case
            1.0 1.0 1.0 ]

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

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

w_shock0 = [ 0.0 0.0 0.0 ;   # no weather shocks - base case
            0.0 0.0 0.0 ]
w = Dict()
for i in 1:length(regions), j in 1:length(states) # weather shock table
    w[regions[i], states[j]] = w_shock0[i,j]
end
################### The Model ##########################

m = MCPModel()               # declare model

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

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

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

#################### Equations #############################
@NLexpression(m, PC[r in regions, s in states],
            (θ * (PE[r,s] * λ[r,s])^(1-esub) + 1-θ)^(1/(1-esub)))
@NLexpression(m, PEU[r in regions],
            sum(πs[s] * PC[r,s]^(1-σ) for s in states)^(1/(1-σ)))
@NLexpression(m, EU[r in regions],
            1/PEU[r])
@NLexpression(m, C[r in regions, s in states],
            EU[r] * (PEU[r]/PC[r,s])^σ)
@NLexpression(m, CM[r in regions],
            K[r] * α[r]/(1-γ) * (1-(1-M[r])^(1-γ)))
#################################################
@mapping(m, 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(m, profit_M[r in regions],
            α[r]/(1-M[r])^γ - sum(πs[s] * PE[r,s] * w[r,s] for s in states))
@mapping(m, profit_X[r in regions, s in states],
            PE[r,s] - PX[s])
@mapping(m, 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(m, market_PR[r in regions],
            0.5 * (1 - K[r]/sqrt(PR[r])))
@mapping(m, market_PX[s in states],
            sum(X[r,s] for r in regions))

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

############## Starting Variable Values ######################
set_start_value.(K, [1.0 1.0])                        #<============== THIS IS WHERE THE ERROR IS
set_start_value.(PE, [1.0 1.0 1.0; 1.0 1.0 1.0])
set_start_value.(PR, [1.0 1.0])
fix.(M, [0.0 0.0]; force = true)
fix.(X, [0.0 0.0 0.0; 0.0 0.0 0.0]; force = true)
fix.(PX, [1.0 1.0 1.0]; force = true)

############## Starting Expression Values ###################
set_start_value.(PEU, [1.0 1.0])                           #<============== CAN I DO THIS
set_start_value.(PC, [1.0 1.0 1.0; 1.0 1.0 1.0])
set_start_value.(EU, [1.0 1.0])
set_start_value.(C, [1.0 1.0 1.0; 1.0 1.0 1.0])
set_start_value.(CM, [1.0 1.0])

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

Here is the Error:

ERROR: LoadError: MethodError: no method matching UnitRange{Int64}(::Array{String,1})
Closest candidates are:
  UnitRange{Int64}(::Any, ::ArrayInterface.StaticInt) where T<:Real at /Users/G/.julia/packages/ArrayInterface/gMtB5/src/static.jl:139
  UnitRange{Int64}(::Any, ::Any) where T<:Real at range.jl:280
  UnitRange{Int64}(::UnitRange{T}) where T<:Real at range.jl:918
chkwon commented 3 years ago

This is related to https://github.com/jump-dev/JuMP.jl/issues/2267.

At this moment, you should use a loop.

njgallagher commented 3 years ago

Thank you for this.