add Multiperiod OPS problem (#15)

Project.toml

@@ -6,18 +6,16 @@ version = "0.2.0"
 [deps]
 InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
-MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Memento = "f28f55f0-a522-5efc-85c2-fe41dfb9b2d9"
 PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
 PowerModelsRestoration = "23adbb12-a187-11e9-26a2-eb4d4e6e68fb"
 
 [compat]
 InfrastructureModels = "~0.6, ~0.7"
-JuMP = "~0.19.1, ~0.20, ~0.21, ~1"
-MathOptInterface = "~0.8, ~0.9, ~1"
-Memento = "~0.10, ~0.11, ~0.12, ~0.13, ~1.0, ~1.1, ~1"
-PowerModels = "~0.18, ~0.19"
-PowerModelsRestoration = "~0.6"
+JuMP = "~0.22, ~0.23, ~1"
+Memento = "~1.0, ~1.1, ~1.2, ~1.3, ~1.4"
+PowerModels = "~0.19.2"
+PowerModelsRestoration = "~0.6, ~0.7"
 julia = "^1.0"
 
 [extras]

src/PowerModelsWildfire.jl

@@ -5,16 +5,15 @@ import InfrastructureModels: nw_id_default
 import PowerModels
 import PowerModelsRestoration
 import JuMP
-import MathOptInterface
 import Memento
 
-const _MOI = MathOptInterface
 const _PM = PowerModels
 const _PMR = PowerModelsRestoration
 const _IM = InfrastructureModels
 
 include("core/variable.jl")
 include("core/constraint_template.jl")
+include("core/constraint.jl")
 include("core/data.jl")
 
 include("form/dcp.jl")

src/core/constraint.jl

@@ -0,0 +1,34 @@
+
+
+""
+function constraint_restoration_indicator_initial(pm::_PM.AbstractPowerModel, n::Int, i::Int)
+    branch = _PM.ref(pm, n, :branch, i)
+    z_branch = _PM.var(pm, n, :z_branch, i)
+    branch_restoration = z_branch = _PM.var(pm, n, :branch_restoration, i)
+
+    JuMP.@constraint(pm.model, branch_restoration <= (1-branch["br_status"]))
+    JuMP.@constraint(pm.model, branch_restoration <= z_branch )
+    JuMP.@constraint(pm.model, branch_restoration >= (1-branch["br_status"]) + z_branch  - 1 )
+end
+
+""
+function constraint_restoration_indicator(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int)
+
+    z_branch_1 = _PM.var(pm, n_1, :z_branch, i)
+    z_branch_2 = _PM.var(pm, n_2, :z_branch, i)
+    branch_restoration =  _PM.var(pm, n_2, :branch_restoration, i)
+
+    JuMP.@constraint(pm.model, branch_restoration <= (1-z_branch_1))
+    JuMP.@constraint(pm.model, branch_restoration <= z_branch_2 )
+    JuMP.@constraint(pm.model, branch_restoration >= (1-z_branch_1) + z_branch_2  - 1 )
+end
+
+""
+function constraint_restoration_budget(pm::_PM.AbstractPowerModel, n::Int, branch_restoration_cost, restoration_budget)
+    branch_restoration =  _PM.var(pm, n, :branch_restoration)
+
+    JuMP.@constraint(pm.model,
+        sum(branch_restoration[id]*cost for (id,cost) in branch_restoration_cost) <= restoration_budget
+    )
+
+end

src/core/constraint_template.jl

@@ -39,3 +39,23 @@ function constraint_branch_active(pm::_PM.AbstractPowerModel, i::Int; nw::Int=nw
     _PMR.constraint_branch_bus_connection(pm, nw, i, branch["t_bus"])
 end
 
+
+""
+function constraint_restoration_budget(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default)
+    if haskey(_PM.ref(pm, nw), :restoration_budget)
+        restoration_budget = _PM.ref(pm, nw, :restoration_budget)
+    else
+        Memento.warn(_PM._LOGGER, "network data should specify restoration_budget, using 10.0 as a default")
+        restoration_budget=10.0
+    end
+    for (id,branch) in _PM.ref(pm, nw, :branch)
+        if !haskey(branch, "restoration_cost")
+            Memento.warn(_PM._LOGGER, "branch data should specify `restoration_cost``, using 10.0 as a default")
+            branch["restoration_cost"]=10.0
+        end
+    end
+
+    branch_restoration_cost = Dict(branch["index"] => branch["restoration_cost"] for (id,branch) in _PM.ref(pm, nw, :branch))
+    constraint_restoration_budget(pm, nw, branch_restoration_cost, restoration_budget)
+end
+

src/core/data.jl

@@ -26,6 +26,19 @@ function calc_load_risk(data)
 end
 
 function calc_risk(data, comp_type)
+    if get(data,"multinetwork",false)==true
+        risk = 0.0
+        for (nwid,nw) in data["nw"]
+            risk += _calc_risk(nw,comp_type)
+        end
+    else
+        risk = _calc_risk(data, comp_type)
+    end
+    return risk
+end
+
+
+function _calc_risk(data, comp_type)
     risk = 0.0
     comp_status = _PM.pm_component_status[comp_type]
     for (comp_id, comp) in get(data,comp_type, Dict())
@@ -41,6 +54,18 @@ function calc_risk(data, comp_type)
 end
 
 function calc_load(data)
+    if get(data,"multinetwork",false)==true
+        load=0.0
+        for (nwid,nw) in data["nw"]
+            load += _calc_load(nw)
+        end
+    else
+        load = _calc_load(data)
+    end
+    return load
+end
+
+function _calc_load(data)
     return sum(load["pd"]*load["status"] for (load_id, load) in data["load"])
 end
 

src/core/export.jl

@@ -3,7 +3,7 @@
 # symbols in your environment, then use `import PowerModelsWildfire` instead of
 # `using PowerModelsWildfire`.
 
-# Do not add PowerModelsRestoration-defined symbols to this exclude list. Instead, rename
+# Do not add PowerModelsWildfire-defined symbols to this exclude list. Instead, rename
 # them with an underscore.
 
 const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]
@@ -24,19 +24,19 @@ end
 # the follow items are also exported for user-friendlyness when calling
 # `using PowerModelsWildfire`
 
-# so that users do not need to import JuMP to use a solver with PowerModels
+# so that users do not need to import JuMP to use a solver with PowerModelsWildfire
 import JuMP: with_optimizer
 export with_optimizer
 
-import MathOptInterface: TerminationStatusCode
+import JuMP: TerminationStatusCode
 export TerminationStatusCode
 
-import MathOptInterface: ResultStatusCode
+import JuMP: ResultStatusCode
 export ResultStatusCode
 
 for status_code_enum in [TerminationStatusCode, ResultStatusCode]
     for status_code in instances(status_code_enum)
-        @eval import MathOptInterface: $(Symbol(status_code))
+        @eval import JuMP: $(Symbol(status_code))
         @eval export $(Symbol(status_code))
     end
 end

src/core/variable.jl

@@ -20,3 +20,50 @@ function variable_bus_active_indicator(pm::_PM.AbstractPowerModel; nw::Int=nw_id
 
     report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :bus, :status, _PM.ids(pm, nw, :bus), z_bus)
 end
+
+
+"variable: `0 <= branch_restoration[l] <= 1` for `l` in `branch`es"
+function variable_branch_restoration_indicator(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true)
+    if relax == false
+        branch_restoration = _PM.var(pm, nw)[:branch_restoration] = JuMP.@variable(pm.model,
+            [l in _PM.ids(pm, nw, :branch)],
+            base_name="$(nw)_branch_restoration",
+            binary = true,
+            start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "branch_restoration_state")
+        )
+    else
+        branch_restoration = _PM.var(pm, nw)[:branch_restoration] = JuMP.@variable(pm.model,
+            [l in _PM.ids(pm, nw, :branch)],
+            base_name="$(nw)_branch_restoration",
+            lower_bound = 0,
+            upper_bound = 1,
+            start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "branch_restoration_state")
+        )
+    end
+    # _PM.var(pm, nw)[:branch_restoration] = z_branch_resto
+
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :branch, :branch_restoration, _PM.ids(pm, nw, :branch), branch_restoration)
+end
+
+"variable: `0 <= load_restoration[l] <= 1` for `l` in `load`es"
+function variable_load_restoration_indicator(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true)
+    if relax == false
+        load_restoration = _PM.var(pm, nw)[:load_restoration] = JuMP.@variable(pm.model,
+            [l in _PM.ids(pm, nw, :load)],
+            base_name="$(nw)_load_restoration",
+            binary = true,
+            start = _PM.comp_start_value(_PM.ref(pm, nw, :load, l), "load_restoration_state")
+        )
+    else
+        load_restoration = _PM.var(pm, nw)[:load_restoration] = JuMP.@variable(pm.model,
+            [l in _PM.ids(pm, nw, :load)],
+            base_name="$(nw)_load_restoration",
+            lower_bound = 0,
+            upper_bound = 1,
+            start = _PM.comp_start_value(_PM.ref(pm, nw, :load, l), "load_restoration_state")
+        )
+    end
+    # _PM.var(pm, nw)[:load_restoration] = load_restoration
+
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :load, :load_restoration, _PM.ids(pm, nw, :load), load_restoration)
+end

src/prob/ops.jl

@@ -1,7 +1,7 @@
 
 ""
 function run_ops(file, model_constructor, optimizer; kwargs...)
-    return _PM.run_model(file, model_constructor, optimizer, build_ops;
+    return _PM.solve_model(file, model_constructor, optimizer, build_ops;
         ref_extensions=[_PM.ref_add_on_off_va_bounds!], kwargs...)
 end
 
@@ -111,3 +111,146 @@ function build_ops(pm::_PM.AbstractPowerModel)
     )
 
 end
+
+""
+function run_mops(file, model_constructor, optimizer, kwargs...)
+    return _PM.solve_model(file, model_constructor, optimizer, build_mn_ops;
+        multinetwork=true, ref_extensions=[_PM.ref_add_on_off_va_bounds!], kwargs...)
+end
+
+
+function build_mn_ops(pm::_PM.AbstractPowerModel)
+    for (n, network) in _PM.nws(pm)
+        variable_bus_active_indicator(pm, nw=n)
+        variable_branch_restoration_indicator(pm, nw=n)
+        variable_load_restoration_indicator(pm, nw=n)
+        _PMR.variable_bus_voltage_on_off(pm, nw=n)
+
+        _PM.variable_gen_indicator(pm, nw=n)
+        _PM.variable_gen_power_on_off(pm, nw=n)
+
+        _PM.variable_storage_indicator(pm, nw=n)
+        _PM.variable_storage_power_mi_on_off(pm, nw=n)
+
+        _PM.variable_branch_indicator(pm, nw=n)
+        _PM.variable_branch_power(pm, nw=n)
+
+        _PM.variable_dcline_power(pm, nw=n)
+
+        _PM.variable_load_power_factor(pm, nw=n, relax=true)
+        _PM.variable_shunt_admittance_factor(pm, nw=n, relax=true)
+
+        _PMR.constraint_model_voltage_damage(pm, nw=n)
+        for i in _PM.ids(pm, :ref_buses, nw=n)
+            _PM.constraint_theta_ref(pm, i, nw=n)
+        end
+
+        for i in _PM.ids(pm, :gen, nw=n)
+            constraint_generation_active(pm, i, nw=n)
+            _PM.constraint_gen_power_on_off(pm, i, nw=n)
+        end
+
+        for i in _PM.ids(pm, :bus, nw=n)
+            constraint_bus_active(pm, i, nw=n)
+            _PMR.constraint_power_balance_shed(pm, i, nw=n)
+        end
+
+        for i in _PM.ids(pm, :storage, nw=n)
+            constraint_storage_active(pm, i, nw=n)
+            _PM.constraint_storage_state(pm, i, nw=n)
+            _PM.constraint_storage_complementarity_mi(pm, i, nw=n)
+            _PM.constraint_storage_on_off(pm,i, nw=n)
+            _PM.constraint_storage_loss(pm, i, nw=n)
+            _PM.constraint_storage_thermal_limit(pm, i, nw=n)
+        end
+
+        for i in _PM.ids(pm, :branch, nw=n)
+            constraint_branch_active(pm, i, nw=n)
+            _PM.constraint_ohms_yt_from_on_off(pm, i, nw=n)
+            _PM.constraint_ohms_yt_to_on_off(pm, i, nw=n)
+
+            _PM.constraint_voltage_angle_difference_on_off(pm, i, nw=n)
+
+            _PM.constraint_thermal_limit_from_on_off(pm, i, nw=n)
+            _PM.constraint_thermal_limit_to_on_off(pm, i, nw=n)
+        end
+
+        for i in _PM.ids(pm, :load, nw=n)
+            constraint_load_active(pm, i, nw=n)
+        end
+
+        for i in _PM.ids(pm, :dcline, nw=n)
+            _PM.constraint_dcline_power_losses(pm, i, nw=n) #not active decision variables
+        end
+
+        constraint_restoration_budget(pm, nw=n)
+    end
+
+    network_ids = sort(collect(_PM.nw_ids(pm)))
+    n_1 = network_ids[1]
+    for i in _PM.ids(pm, :branch, nw=n_1)
+        constraint_restoration_indicator_initial(pm, n_1, i)
+    end
+
+    for n_2 in network_ids[2:end]
+        for i in _PM.ids(pm, :branch, nw=n_2)
+            constraint_restoration_indicator(pm, n_1, n_2, i)
+        end
+        n_1 = n_2
+    end
+
+    # Add Objective Function
+    # ----------------------
+    # Maximize power delivery while minimizing wildfire risk
+    n_1 = network_ids[1]
+    if haskey(_PM.ref(pm, n_1), :risk_weight)
+        alpha = _PM.ref(pm, n_1, :risk_weight)
+    else
+        Memento.warn(_PM._LOGGER, "network data should specify risk_weight, using 0.5 as a default")
+        alpha = 0.5
+    end
+
+    if haskey(_PM.ref(pm, n_1), :disable_cost)
+        disable_cost = _PM.ref(pm, n_1, :disable_cost)
+    else
+        Memento.warn(_PM._LOGGER, "network data should specify disable_cost, using 10.0 as a default")
+        disable_cost = 10.0
+    end
+
+
+
+    for comp_type in [:gen, :load, :bus, :branch]
+        for nwid in _PM.nw_ids(pm)
+            for (id,comp) in  _PM.ref(pm, nwid, comp_type)
+                if ~haskey(comp, "power_risk")
+                    @warn "$(comp_type) $(id) does not have a power_risk value, using 0.0 as a default"
+                    comp["power_risk"] = 0.0
+                end
+            end
+        end
+    end
+
+    # scale based on total load demand and risk
+    total_load = sum(sum(load["pd"] for (load_id,load) in nw[:load]) for (nwid,nw) in  _PM.nws(pm))
+    total_risk =
+    sum(sum(sum(get(comp,"power_risk",0)
+                for (compid,comp) in  _PM.ref(pm, nwid, comp_type))
+            for nwid in  _PM.nw_ids(pm))
+        for comp_type in [:branch,:gen,:bus,:load]
+    )
+
+    z_demand = Dict(nwid => _PM.var(pm, nwid, :z_demand) for nwid in _PM.nw_ids(pm))
+    z_branch = Dict(nwid => _PM.var(pm, nwid, :z_branch) for nwid in _PM.nw_ids(pm))
+
+    JuMP.@objective(pm.model, Max,
+        sum(
+            (1-alpha)*(
+                sum(load["pd"]*z_demand[nwid][i]/total_load for (i,load) in _PM.ref(pm, nwid, :load))
+            )
+            -alpha*(
+                sum(z_branch[nwid][i]*branch["power_risk"]/total_risk for (i,branch) in _PM.ref(pm, nwid, :branch))
+                +sum((1-z_branch[nwid][i])*disable_cost/total_risk for (i,branch) in _PM.ref(pm, nwid, :branch))
+            )
+        for nwid in _PM.nw_ids(pm))
+    )
+end