Skip to content

Instantly share code, notes, and snippets.

@pierre-haessig

pierre-haessig/JuMP_modular_problem.ipynb

Last active July 30, 2024 08:33
Show Gist options
  • Save pierre-haessig/79358e7c45e7ab9954938188162d92b6 to your computer and use it in GitHub Desktop.
Save pierre-haessig/79358e7c45e7ab9954938188162d92b6 to your computer and use it in GitHub Desktop.
Download ZIP
Modular implementation of optimization problems with JuMP
Display the source blob
Display the rendered blob
Raw
JuMP_modular_problem.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "66d1b45f-a4dc-458c-99f2-6cf1f0d0ef05",
"metadata": {},
"source": [
"# Modular implementation of optimization problems with JuMP\n",
"\n",
"An exploration of **modular implementations** for two-stage stochastic optimization of energy systems.\n",
"\n",
"Key question is how to achieve **modularity**, in particular having about the same model definition to handle the following variants of energy system planning:\n",
"1. single stage deterministic variant (1 set of sizing variables + 1 set operation variables)\n",
"2. two stage deterministic variant (two copies of the single stage problem, the first one for the first years of a project and the second one being a “resize” at mid-project)\n",
"3. two-stage stochastic variant ($1+n_s$ copies of the single stage problem, with $n_s$ the number of stochastic scenarios for the second period of the project)\n",
"\n",
"Modularity is welcomed because each stage has more or less the same variables. However, in JuMP, variables declared with the `@variable` macro are unique. Storing variables in an array is acceptable for the second stage variables, but not for first and second stage since the dimensions are not the same. Also, the variant 1 (only one stage) has to be supported.\n",
"\n",
"**Content**:\n",
"\n",
"- This notebook contains a dummy energy system optimization problem to test modular implementations.\n",
"- In the first sections, a classical \"flat\" JuMP formulation is given to specify the expected result.\n",
"- Then, a modular version is provided. In the design process, I ended up with two versions (\"A\" and \"B\") and in the end I prefer option B because it creates less noise in the inner model.\n",
"\n",
"Pierre Haessig, June 2024, [CC-BY 4.0](href=\"https://creativecommons.org/licenses/by/4.0/)\n",
"\n",
"- notebook posted on Github Gist: https://gist.github.com/pierre-haessig/79358e7c45e7ab9954938188162d92b6\n",
"- companion post on Forum https://discourse.julialang.org/t/modular-implementation-of-optimization-problems-with-jump/116261\n",
"\n",
"---\n",
"\n",
"*Late edit* (July 30, 2024): following the latest discussion development on the forum, I tend to prefer not using an customized `add_var` routine and instead use the regular `@variable`.\n",
"\n",
"See fragment from my July 17th post:\n",
"\n",
"```julia\n",
"data[\"x\"] = @variable(model, 0 <= x <= 1)\n",
"data[\"y\"] = @variable(model, 0 <= y <= 1)\n",
"@constraint(model, x+y <= 1) # here comes the interest of having x and y in the namespace\n",
"# ...\n",
"\n",
"for var in keys(object_dictionary(model))\n",
" unregister(model, var)\n",
"end\n",
"```\n",
"\n",
"In the last iteration, the extra burden is quite low:\n",
"- prefix each line with `data[\"varname\"] = `\n",
"- unregister all symbols from the model at the end (3 lines)\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "2c15f168-2fa1-4221-aa40-165e6c9ba9c6",
"metadata": {},
"source": [
"## Docs and references\n",
"\n",
"Discussions e.g. on dynamic generation of JuMP variable names: \n",
"\n",
"- First, background knowledge: JuMP doc clarifying the different notions of \"names\" with JuMP variables: [String names, symbolic names, and bindings](https://jump.dev/JuMP.jl/stable/manual/variables/#variable_names_and_bindings)\n",
"- Discourse “How to merge two `JuMP` models?” https://discourse.julialang.org/t/how-to-merge-two-jump-models/113646, 2024\n",
" - odow points to [Design patterns for larger models](https://jump.dev/JuMP.jl/stable/tutorials/getting_started/design_patterns_for_larger_models/) where problem variables constraints are added within function, but *still the variable names* are harcoded\n",
"- Stack Overflow: Impossibility to dynamically generate JuMP variable names: “in no case can you dynamically create a binding like x_i [with k=1:K]”, 2022 https://stackoverflow.com/questions/72692238/define-several-jump-variables-names-in-loop-based-on-iteration, 2022\n",
"- Discourse: “Using anonymous and or symbols to create dynamic varible names?”\n",
" - “you should create an *appropriate data structure* to hold the variables (`Vector`, a `Dict`, or some user-defined `struct`: the choice is problem-dependent) https://discourse.julialang.org/t/using-anonymous-and-or-symbols-to-create-dynamic-varible-names/83095/2, 2022\n",
" - I posted the example with `add_stage_variable!` below to clarify my interrogations...\n",
" - odow's answer: only the dynamic generation of Julia variables binding is discouraged. Dynamic generation of JuMP registered variables is fine \n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "377d6d0c-2be5-4012-b8a5-7228962d3ecc",
"metadata": {},
"outputs": [],
"source": [
"using JuMP\n",
"using HiGHS"
]
},
{
"cell_type": "markdown",
"id": "442e8764-48fa-4961-9af9-b1730198884f",
"metadata": {},
"source": [
"## Experiment: a dynamic `add_stage_variable!` function\n",
"\n",
"Posted on Discourse for feedback: https://discourse.julialang.org/t/using-anonymous-and-or-symbols-to-create-dynamic-varible-names/83095/3. Positive feedback received from Oscar Dowson.\n",
"\n",
"Remark:\n",
"`data1` gets printed as `Dict{String, Any}(\"x\" => x1)`. However, by modifying `add_stage_variable!`, it could have been `\"x1\" => x1`, `\"x1\" => x` or `\"x\" => x` depending on taste, since there is no risk of clash with the content of `data2`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "20511c32-bccd-4328-afa1-0d2e41a9936f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Dict{String, Any}(\"x\" => x1)\n"
]
},
{
"data": {
"text/plain": [
"A JuMP Model\n",
"Feasibility problem with:\n",
"Variables: 2\n",
"Model mode: AUTOMATIC\n",
"CachingOptimizer state: EMPTY_OPTIMIZER\n",
"Solver name: HiGHS\n",
"Names registered in the model: x1, x2"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"add variable `name` to model `m` and to Dict `dat`, for optional `stage`\"\n",
"function add_stage_variable_test!(m::Model, dat::Dict{String,Any}, name::String; stage=\"\")\n",
" name_suffix = \"$name$stage\"\n",
" dat[name] = @variable(m, base_name=name_suffix) # here I can use name or name_suffix\n",
" # Optional last line: finish the job by registering variable as name_suffix????\n",
" m[Symbol(name_suffix)] = dat[name]\n",
"end\n",
"\n",
"m = Model(HiGHS.Optimizer)\n",
"\n",
"data1 = Dict{String,Any}()\n",
"data2 = Dict{String,Any}()\n",
"add_stage_variable_test!(m, data1, \"x\"; stage=1)\n",
"add_stage_variable_test!(m, data2, \"x\"; stage=2)\n",
"\n",
"println(data1)\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "84a03daf-5a9b-4e5e-8221-03c781c6e89c",
"metadata": {},
"source": [
"### About the difference between String names and registered names\n",
"\n",
"It's allowed to have two different (anonymous or not) variable with the *same name attribute*. They are just *not* the same variable...\n",
"See [String names, symbolic names, and bindings](https://jump.dev/JuMP.jl/stable/manual/variables/#variable_names_and_bindings)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "75258298-2376-48be-9139-5ea0fc52e2c7",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(x + x, 2 x, 2 x)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = Model(HiGHS.Optimizer)\n",
"\n",
"@variable(m, x)\n",
"x2 = @variable(m, base_name=\"x\")\n",
"\n",
"x+x2, x+x, x2+x2 # x!=x2 despite their string name, which can be confusing!"
]
},
{
"cell_type": "markdown",
"id": "7fdb63da-a8ed-42d6-948d-d286618621ca",
"metadata": {},
"source": [
"However, adding variable with same to-be-registered name is not possible:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "21abe4b2-09bd-4f47-b780-52c79306b33b",
"metadata": {},
"outputs": [
{
"ename": "LoadError",
"evalue": "An object of name x is already attached to this model. If this\n is intended, consider using the anonymous construction syntax, for example,\n `x = @variable(model, [1:N], ...)` where the name of the object does\n not appear inside the macro.\n\n Alternatively, use `unregister(model, :x)` to first unregister\n the existing name from the model. Note that this will not delete the\n object; it will just remove the reference at `model[:x]`.\n",
"output_type": "error",
"traceback": [
"An object of name x is already attached to this model. If this\n is intended, consider using the anonymous construction syntax, for example,\n `x = @variable(model, [1:N], ...)` where the name of the object does\n not appear inside the macro.\n\n Alternatively, use `unregister(model, :x)` to first unregister\n the existing name from the model. Note that this will not delete the\n object; it will just remove the reference at `model[:x]`.\n",
"",
"Stacktrace:",
" [1] error(s::String)",
" @ Base ./error.jl:35",
" [2] _error_if_cannot_register(model::Model, name::Symbol)",
" @ JuMP ~/.julia/packages/JuMP/7rBNn/src/macros.jl:411",
" [3] macro expansion",
" @ ~/.julia/packages/JuMP/7rBNn/src/macros.jl:400 [inlined]",
" [4] top-level scope",
" @ In[4]:1"
]
}
],
"source": [
"@variable(m, x)"
]
},
{
"cell_type": "markdown",
"id": "f5933cd8-1dd9-4abe-85db-5d415c9c974f",
"metadata": {},
"source": [
"Registered variables (and constraints...) are accessible with `object_dictionary`:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "746fc76a-660b-4637-89e6-e9a8eee0fbdb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{Symbol, Any} with 1 entry:\n",
" :x => x"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"object_dictionary(m)"
]
},
{
"cell_type": "markdown",
"id": "5cae30db-6947-4ce0-ac8c-6565ca95c79c",
"metadata": {},
"source": [
"## Optimization problem definition using flat implementations\n",
"\n",
"### Warm-up: a single stage problem\n",
"\n",
"A dummy energy system optimization model to explore the different implementations: only one power source \"gen\" to feed a load on `K` time instant. The generator has rated capacity `Pgen` to be optimized. The solution is known in advance: Pgen = max(load)."
]
},
{
"cell_type": "markdown",
"id": "eb5248d0-22f1-44b1-9c34-a9b09edf3612",
"metadata": {},
"source": [
"Problem data: time series of the load to be supplied along time"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "dfad88f6-3b5b-4d1a-9bb3-8cb2d377298a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time length: 4\n"
]
}
],
"source": [
"load = [1., 1.5, 1.2, 0.9]\n",
"K = length(load)\n",
"println(\"time length: $K\")"
]
},
{
"cell_type": "markdown",
"id": "488dfa6f-7b9a-4124-bec0-357442f3f6ff",
"metadata": {},
"source": [
"#### Flat implementation of single stage"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "20e142a4-e858-47e5-815d-3a570d728b2b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pgen=1.5\n"
]
}
],
"source": [
"m = Model(HiGHS.Optimizer)\n",
"set_silent(m)\n",
"\n",
"@variable(m, Pgen)\n",
"@variable(m, gen[1:K])\n",
"\n",
"@constraint(m, balance[k=1:K], gen[k] == load[k])\n",
"@constraint(m, gen .<= Pgen)\n",
"\n",
"cost = Pgen + sum(gen)\n",
"@objective(m, Min, cost)\n",
"\n",
"optimize!(m)\n",
"println(\"Pgen=$(value(Pgen))\")"
]
},
{
"cell_type": "markdown",
"id": "1882e358-9c75-4b1d-a16e-09983c98123b",
"metadata": {},
"source": [
"#### Single stage wrapped-in-a-function version\n",
"\n",
"The function allows reuse, but since the variable names are unique, it cannot be used to instantiate multi stages/scenarios"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "dfea52dc-ba99-4aac-a0cf-55129b4e085a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pgen=3.0\n"
]
}
],
"source": [
"function optim_model!(m, load)\n",
" K = length(load)\n",
" @variable(m, Pgen)\n",
" @variable(m, gen[1:K])\n",
" \n",
" @constraint(m, balance[k=1:K], gen[k] == load[k])\n",
" @constraint(m, gen .<= Pgen)\n",
" \n",
" @objective(m, Min, Pgen+sum(gen))\n",
"end\n",
"\n",
"m = Model(HiGHS.Optimizer)\n",
"set_silent(m)\n",
"\n",
"# Optim with twice the load to see the difference:\n",
"optim_model!(m, load*2)\n",
"\n",
"optimize!(m)\n",
"println(\"Pgen=$(value(m[:Pgen]))\")"
]
},
{
"cell_type": "markdown",
"id": "ba265314-75ca-42ed-99b1-01fb404d89b5",
"metadata": {},
"source": [
"### Two-stage Stochastic Program"
]
},
{
"cell_type": "markdown",
"id": "eb87aecf-a353-49a6-b0f8-2600d75261f9",
"metadata": {},
"source": [
"#### Data for 2nd stage scenarios\n",
"\n",
"3 equally probable load changes"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "43365b27-6f56-4251-8eeb-75f08f441208",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of scenarios: 3\n"
]
}
],
"source": [
"load_scenar = [0.5, 1.0, 2.0] # (decrease, constant, increase)\n",
"ns = length(load_scenar)\n",
"proba_scenar = [1., 1., 1.]/ns\n",
"println(\"number of scenarios: $ns\")"
]
},
{
"cell_type": "markdown",
"id": "a3b4186f-4794-4bde-adb1-d7cf2e51e7a7",
"metadata": {},
"source": [
"#### Flat formulation of the two-stage SP\n",
"\n",
"to specify the expected results of subsequent implementations\n",
"\n",
"To make the problem a bit interesting, there is some coupling between the two stages. Here the rated powers at second stage must be >= than the rated power at first stage."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e6943b1f-05c0-4cb8-8f29-cabd9aa8554b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:\n",
" Pgen1 - Pgen2[1] ≤ 0\n",
" Pgen1 - Pgen2[2] ≤ 0\n",
" Pgen1 - Pgen2[3] ≤ 0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = Model(HiGHS.Optimizer)\n",
"set_silent(m)\n",
"\n",
"@variable(m, Pgen1)\n",
"@variable(m, Pgen2[1:ns])\n",
"\n",
"@variable(m, gen1[1:K])\n",
"@variable(m, gen2[1:K, 1:ns]) # timeseries scenarios stacked as columns\n",
"\n",
"@constraint(m, balance1[k=1:K], gen1[k] == load[k])\n",
"@constraint(m, balance2[k=1:K, s=1:ns], gen2[k,s] == load[k]*load_scenar[s])\n",
"@constraint(m, gen1 .<= Pgen1)\n",
"@constraint(m, [s=1:ns], gen2[:,s] .<= Pgen2[s])\n",
"\n",
"# Coupling between the two stages\n",
"@constraint(m, Pgen1 .<= Pgen2)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "0a3992c1-dc72-4e25-b6b4-008397e719d9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"A JuMP Model\n",
"Minimization problem with:\n",
"Variables: 20\n",
"Objective function type: AffExpr\n",
"`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 16 constraints\n",
"`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 19 constraints\n",
"Model mode: AUTOMATIC\n",
"CachingOptimizer state: EMPTY_OPTIMIZER\n",
"Solver name: HiGHS\n",
"Names registered in the model: Pgen1, Pgen2, balance1, balance2, gen1, gen2"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Costs\n",
"cost1 = Pgen1 + sum(gen1)\n",
"\n",
"Pgen2_exp = Pgen2' * proba_scenar\n",
"gen2_exp = gen2 * proba_scenar # inner product at each instant=line\n",
"cost2 = Pgen2_exp + sum(gen2_exp)\n",
"\n",
"@objective(m, Min, cost1 + cost2)\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "4d77375b-9b63-47a9-864b-8715dc391db7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pgen1=1.5\n",
"Pgen2=[1.5, 1.5, 3.0]\n"
]
}
],
"source": [
"optimize!(m)\n",
"println(\"Pgen1=$(value(Pgen1))\")\n",
"println(\"Pgen2=$(value.(Pgen2))\")"
]
},
{
"cell_type": "markdown",
"id": "cd0ede44-f844-4111-8d59-8b6a7a20ed29",
"metadata": {},
"source": [
"## Modular single stage / two-stage implementations"
]
},
{
"cell_type": "markdown",
"id": "8bd2ad3b-f584-4489-b6a1-7d75729dcbc5",
"metadata": {},
"source": [
"### Modular implementation A\n",
"\n",
"A rather compact implementation of the problem which can be used for a single stage but also for multiple stages.\n",
"In this implementation, the lower level model can be composed of a **single or multiple scenarios**.\n",
"The two stage aspect is implemented in an outer composition of **two inner model instances**.\n",
"Main properties:\n",
"\n",
"- objective is left undefined so that it can be set outside, but the stage cost is made available in `data[\"cost\"]`\n",
"- variable definition is compact thanks to `add_stage_variable!`\n",
"- however, constraint and cost definitions are not compact: there is still a need to switch depending on the number of scenarios (scalar vs vector).\n",
" - for the cost, an aggregation function (Expectation, CVaR...) could be implemented quite easily\n",
" - but for constraints, it is more tricky because several subcases need to be handled: type of constraint (`<=`, `==`, `>=`) and different broadcastings depending on the two sides (times series or scalar...)"
]
},
{
"cell_type": "markdown",
"id": "04c097c6-0a87-4779-80b1-f03417ffd4fb",
"metadata": {},
"source": [
"#### Helpers version A"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "cff8c0d4-848d-4934-989d-1e24689aa986",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"add_stage_variable!"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"add variable `name` to model `m` and to Dict `dat`, for optional `stage`.\n",
" \n",
"The variable is optionnally a Vector of length `K` or `ns` or Array of shape (`K`,`ns`)\n",
"if `K` and/or `ns` >1. (timeseries scenarios stacked as columns)\n",
"\n",
"returns the variable\n",
"\"\"\"\n",
"function add_stage_variable!(m::Model, dat::Dict{String,Any}, name::String;\n",
" K::Int=1, stage=\"\", ns::Int=1)\n",
" name_suffix = \"$name$stage\"\n",
" if ns==1 && K==1\n",
" v = @variable(m, base_name=name_suffix)\n",
" elseif ns==1 && K>1\n",
" v = @variable(m, [1:K], base_name=name_suffix)\n",
" elseif ns>1 && K==1\n",
" v = @variable(m, [1:ns], base_name=name_suffix)\n",
" elseif ns>1 && K>1\n",
" v = @variable(m, [1:K, 1:ns], base_name=name_suffix)\n",
" else\n",
" throw(ArgumentError(\"`ns` and `K` array size parameters should be >=1\"))\n",
" end\n",
" dat[name] = v\n",
" return v\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "7e4a02ce-0d2f-45f2-ab6a-8b69823cc6cd",
"metadata": {},
"source": [
"#### Inner model A (single stage with one or several scenarios)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "6d05ca3c-fc67-4a95-8a10-050aa009bce1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"optim_model_stage! (generic function with 1 method)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function optim_model_stage!(m, dat, load; stage=\"\", load_scenar=[1.0], proba_scenar=[1.0])\n",
" K = length(load)\n",
" ns = length(load_scenar)\n",
" @assert length(load_scenar) == length(proba_scenar)\n",
"\n",
" # Variables\n",
" Pgen = add_stage_variable!(m, dat, \"Pgen\"; stage=stage, ns=ns)\n",
" gen = add_stage_variable!(m, dat, \"gen\"; K=K, stage=stage, ns=ns) # timeseries scenarios stacked as columns\n",
" \n",
" # Balance named constraint: TODO create a add_stage_constraint function which handles scenarios?\n",
" dat[\"balance\"] = ns==1 ?\n",
" @constraint(m, [k=1:K], gen[k] == load[k], base_name=\"balance$stage\") :\n",
" @constraint(m, [k=1:K, s=1:ns], gen[k,s] == load[k]*load_scenar[s], base_name=\"balance$stage\")\n",
"\n",
" dat[\"Pgen_cons\"] = ns==1 ?\n",
" @constraint(m, [k=1:K], gen[k] <= Pgen) :\n",
" @constraint(m, [k=1:K, s=1:ns], gen[k,s] <= Pgen[s])\n",
"\n",
" dat[\"cost\"] = if ns == 1\n",
" Pgen + sum(gen)\n",
" else\n",
" # alternative: Pgen_exp = add_stage_variable!(m, dat, \"Pgen_exp\"; stage=stage)\n",
" # constraint Pgen_exp == Pgen' * proba_scenar\n",
" dat[\"Pgen_exp\"] = Pgen' * proba_scenar\n",
" dat[\"gen_exp\"] = gen * proba_scenar # inner product at each instant=line\n",
" dat[\"Pgen_exp\"] + sum(dat[\"gen_exp\"])\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "b0ff75b2-5c3b-47cc-82d1-b394b22d752d",
"metadata": {},
"source": [
"#### Single stage usage"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "c6ce3f3d-fce6-4405-b807-d0b998a31f9b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{String, Any} with 5 entries:\n",
" \"balance\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"cost\" => Pgen + gen[1] + gen[2] + gen[3] + gen[4]\n",
" \"gen\" => VariableRef[gen[1], gen[2], gen[3], gen[4]]\n",
" \"Pgen_cons\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"Pgen\" => Pgen"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = Model(HiGHS.Optimizer)\n",
"set_silent(m)\n",
"\n",
"data = Dict{String,Any}()\n",
"\n",
"optim_model_stage!(m, data, load)\n",
"\n",
"data"
]
},
{
"cell_type": "markdown",
"id": "a842c4a1-a18e-4325-8649-5124d2c1d06c",
"metadata": {},
"source": [
"Set optimization objective"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "b6b90cc4-29d1-4687-bf0b-9f5d77133496",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"A JuMP Model\n",
"Minimization problem with:\n",
"Variables: 5\n",
"Objective function type: AffExpr\n",
"`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 4 constraints\n",
"`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 4 constraints\n",
"Model mode: AUTOMATIC\n",
"CachingOptimizer state: EMPTY_OPTIMIZER\n",
"Solver name: HiGHS"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@objective(m, Min, data[\"cost\"])\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "4c205822-343f-4dc8-aa9c-1656b731c692",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pgen=1.5\n"
]
}
],
"source": [
"optimize!(m)\n",
"Pgen_val = value(data[\"Pgen\"])\n",
"println(\"Pgen=$(Pgen_val)\")"
]
},
{
"cell_type": "markdown",
"id": "9ab6bda4-b21e-4288-867d-db387a2d49ae",
"metadata": {},
"source": [
"#### Two-stage usage\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "91c8d7b3-72c8-413b-8a91-2ab7ff3c4015",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{String, Any} with 7 entries:\n",
" \"balance\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"cost\" => 0.3333333333333333 Pgen2[1] + 0.3333333333333333 Pgen2[2] + 0.…\n",
" \"gen\" => VariableRef[gen2[1,1] gen2[1,2] gen2[1,3]; gen2[2,1] gen2[2,2]…\n",
" \"Pgen_cons\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"gen_exp\" => AffExpr[0.3333333333333333 gen2[1,1] + 0.3333333333333333 gen2…\n",
" \"Pgen_exp\" => 0.3333333333333333 Pgen2[1] + 0.3333333333333333 Pgen2[2] + 0.…\n",
" \"Pgen\" => VariableRef[Pgen2[1], Pgen2[2], Pgen2[3]]"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = Model(HiGHS.Optimizer)\n",
"set_silent(m)\n",
"\n",
"data1 = Dict{String,Any}()\n",
"data2 = Dict{String,Any}()\n",
"\n",
"optim_model_stage!(m, data1, load; stage=1)\n",
"optim_model_stage!(m, data2, load; stage=2,\n",
" load_scenar=load_scenar, proba_scenar=proba_scenar)\n",
"\n",
"data2"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "725256a0-ebdd-467b-9ac5-2537a72c836b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$ 0.3333333333333333 Pgen2_{1} + 0.3333333333333333 Pgen2_{2} + 0.3333333333333333 Pgen2_{3} + 0.3333333333333333 gen2_{1,1} + 0.3333333333333333 gen2_{1,2} + 0.3333333333333333 gen2_{1,3} + 0.3333333333333333 gen2_{2,1} + 0.3333333333333333 gen2_{2,2} + 0.3333333333333333 gen2_{2,3} + 0.3333333333333333 gen2_{3,1} + 0.3333333333333333 gen2_{3,2} + 0.3333333333333333 gen2_{3,3} + 0.3333333333333333 gen2_{4,1} + 0.3333333333333333 gen2_{4,2} + 0.3333333333333333 gen2_{4,3} $"
],
"text/plain": [
"0.3333333333333333 Pgen2[1] + 0.3333333333333333 Pgen2[2] + 0.3333333333333333 Pgen2[3] + 0.3333333333333333 gen2[1,1] + 0.3333333333333333 gen2[1,2] + 0.3333333333333333 gen2[1,3] + 0.3333333333333333 gen2[2,1] + 0.3333333333333333 gen2[2,2] + 0.3333333333333333 gen2[2,3] + 0.3333333333333333 gen2[3,1] + 0.3333333333333333 gen2[3,2] + 0.3333333333333333 gen2[3,3] + 0.3333333333333333 gen2[4,1] + 0.3333333333333333 gen2[4,2] + 0.3333333333333333 gen2[4,3]"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data2[\"cost\"]"
]
},
{
"cell_type": "markdown",
"id": "1bab0e3e-c01a-4bed-8197-b18fc5cdb1b4",
"metadata": {},
"source": [
"Coupling between the two stages"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "affe4401-5b20-4ff1-98b3-4cd9728440fc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:\n",
" Pgen1 - Pgen2[1] ≤ 0\n",
" Pgen1 - Pgen2[2] ≤ 0\n",
" Pgen1 - Pgen2[3] ≤ 0"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@constraint(m, data1[\"Pgen\"] .<= data2[\"Pgen\"])"
]
},
{
"cell_type": "markdown",
"id": "8fd45689-aff7-4c31-9870-462eec98884b",
"metadata": {},
"source": [
"Set optimization objective, summing both stages"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "eefdc82f-16a7-47d1-8f36-08e3447d22bf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"A JuMP Model\n",
"Minimization problem with:\n",
"Variables: 20\n",
"Objective function type: AffExpr\n",
"`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 16 constraints\n",
"`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 19 constraints\n",
"Model mode: AUTOMATIC\n",
"CachingOptimizer state: EMPTY_OPTIMIZER\n",
"Solver name: HiGHS"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@objective(m, Min, data1[\"cost\"] + data2[\"cost\"])\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "7919f901-0cd4-4d77-9776-d1fafc495985",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pgen1=1.5\n",
"Pgen2=[1.5, 1.5, 3.0]\n"
]
}
],
"source": [
"optimize!(m)\n",
"Pgen1_val = value(data1[\"Pgen\"])\n",
"println(\"Pgen1=$(Pgen1_val)\")\n",
"Pgen2_val = value.(data2[\"Pgen\"])\n",
"println(\"Pgen2=$(Pgen2_val)\")"
]
},
{
"cell_type": "markdown",
"id": "68d0c212-1850-4e26-885f-7bb7da387251",
"metadata": {},
"source": [
"### Modular implementation B\n",
"\n",
"In this implementation, the lower level model implements only a **single scenario**, which makes for a much lighter inner model compared to version A.\n",
"The two stage aspect is implemented in an outer composition of $(1+ns)$ inner model instances, so there is more work left of the outer model.\n",
"\n",
"In particular, computing the expected cost is left for the outer model.\n",
"\n",
"Main properties:\n",
"\n",
"- objective is left undefined so that it can be set outside, but the scenario cost is made available in `data[\"cost\"]`"
]
},
{
"cell_type": "markdown",
"id": "a612fbc1-04ca-43bf-8398-e4ee9f8e2bd8",
"metadata": {},
"source": [
"#### Helpers version B"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "b12d67fb-4601-4ace-b787-23bf0bc8d83f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"add_stage_variable_B!"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"add variable `name` to model `m` and to Dict `dat`, for optional `stage`.\n",
" \n",
"The variable is optionnally a Vector of length `K` if `K` >1 (timeseries).\n",
"\n",
"`stage` is used as a suffix for the `base_name` of the JuMP variable\n",
"\n",
"returns the variable\n",
"\"\"\"\n",
"function add_stage_variable_B!(m::Model, dat::Dict{String,Any}, name::String;\n",
" K::Int=1, stage=\"\")\n",
" name_suffix = \"$name$stage\"\n",
" if K==1\n",
" v = @variable(m, base_name=name_suffix)\n",
" elseif K>1\n",
" v = @variable(m, [1:K], base_name=name_suffix)\n",
" else\n",
" throw(ArgumentError(\"`K` array size parameter should be >=1\"))\n",
" end\n",
" dat[name] = v\n",
" return v\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "789b1c1e-ee52-499f-816b-e35a11716382",
"metadata": {},
"source": [
"#### Inner model B (single scenario)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "70ef5b0a-243d-4391-b2f2-d41a2fe27fd3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"optim_model_stage_B! (generic function with 1 method)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function optim_model_stage_B!(m, dat, load; stage=\"\")\n",
" K = length(load)\n",
"\n",
" # Variables\n",
" Pgen = add_stage_variable_B!(m, dat, \"Pgen\"; stage=stage)\n",
" gen = add_stage_variable_B!(m, dat, \"gen\"; K=K, stage=stage)\n",
" \n",
" # Constraints:\n",
" dat[\"balance\"] = @constraint(m, [k=1:K], gen[k] == load[k], base_name=\"balance$stage\")\n",
" dat[\"Pgen_cons\"] = @constraint(m, [k=1:K], gen[k] <= Pgen, base_name=\"Pgen_cons$stage\")\n",
"\n",
" # scenario cost\n",
" dat[\"cost\"] = Pgen + sum(gen)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "01a26f96-09ff-425d-b1dc-72378d9df8a1",
"metadata": {},
"source": [
"#### Single stage usage"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "a25e50c4-42a7-41ea-a957-6995f2c48e96",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{String, Any} with 5 entries:\n",
" \"balance\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"cost\" => Pgen + gen[1] + gen[2] + gen[3] + gen[4]\n",
" \"gen\" => VariableRef[gen[1], gen[2], gen[3], gen[4]]\n",
" \"Pgen_cons\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"Pgen\" => Pgen"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = Model(HiGHS.Optimizer)\n",
"set_silent(m)\n",
"\n",
"data = Dict{String,Any}()\n",
"\n",
"optim_model_stage_B!(m, data, load)\n",
"\n",
"data"
]
},
{
"cell_type": "markdown",
"id": "0e75e099-e8f0-4d4c-99ac-0a82628d8f6d",
"metadata": {},
"source": [
"Set optimization objective"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "312201d8-ae31-4b01-a6ec-e65e23faacdb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"A JuMP Model\n",
"Minimization problem with:\n",
"Variables: 5\n",
"Objective function type: AffExpr\n",
"`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 4 constraints\n",
"`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 4 constraints\n",
"Model mode: AUTOMATIC\n",
"CachingOptimizer state: EMPTY_OPTIMIZER\n",
"Solver name: HiGHS"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@objective(m, Min, data[\"cost\"])\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "987b447e-c440-424e-9f30-e0e2cf17ab4a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pgen=1.5\n"
]
}
],
"source": [
"optimize!(m)\n",
"Pgen_val = value(data[\"Pgen\"])\n",
"println(\"Pgen=$(Pgen_val)\")"
]
},
{
"cell_type": "markdown",
"id": "615474b8-f1c1-4700-a612-5776fc968d10",
"metadata": {},
"source": [
"#### Two-stage usage\n"
]
},
{
"cell_type": "markdown",
"id": "1100d6f9-0852-4240-8d41-b03a6a070dc9",
"metadata": {},
"source": [
"Compared to Version A, the `data2` is a Vector of Dict, not just a Dict. Then, the outer model needs to collect the cost of each scenario of the second stage to compute an aggregated (e.g. expected) costs."
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "00634ed5-edfd-46c1-a3b9-3a094c7db375",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{String, Any} with 5 entries:\n",
" \"balance\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"cost\" => Pgen21 + gen21[1] + gen21[2] + gen21[3] + gen21[4]\n",
" \"gen\" => VariableRef[gen21[1], gen21[2], gen21[3], gen21[4]]\n",
" \"Pgen_cons\" => ConstraintRef{Model, ConstraintIndex{ScalarAffineFunction{Floa…\n",
" \"Pgen\" => Pgen21"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = Model(HiGHS.Optimizer)\n",
"set_silent(m)\n",
"\n",
"data1 = Dict{String,Any}()\n",
"optim_model_stage_B!(m, data1, load; stage=1)\n",
"\n",
"data2 = [Dict{String,Any}() for s=1:ns]\n",
"for s=1:ns\n",
" load_s = load*load_scenar[s]\n",
" optim_model_stage_B!(m, data2[s], load_s; stage=\"2$s\")\n",
"end\n",
"\n",
"data2[1]"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "f518732b-5a8a-4345-bb5e-48d71c0eb46f",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$ Pgen21 + gen21_{1} + gen21_{2} + gen21_{3} + gen21_{4} $"
],
"text/plain": [
"Pgen21 + gen21[1] + gen21[2] + gen21[3] + gen21[4]"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data2[1][\"cost\"]"
]
},
{
"cell_type": "markdown",
"id": "198be7da-9632-42ae-b6c9-10f2a9142c89",
"metadata": {},
"source": [
"Creating expected cost of second stage (extra work for outer model compared to option A)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "29dc4e8a-9420-4dd7-9db1-7c3adc52626b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$ 0.3333333333333333 Pgen21 + 0.3333333333333333 gen21_{1} + 0.3333333333333333 gen21_{2} + 0.3333333333333333 gen21_{3} + 0.3333333333333333 gen21_{4} + 0.3333333333333333 Pgen22 + 0.3333333333333333 gen22_{1} + 0.3333333333333333 gen22_{2} + 0.3333333333333333 gen22_{3} + 0.3333333333333333 gen22_{4} + 0.3333333333333333 Pgen23 + 0.3333333333333333 gen23_{1} + 0.3333333333333333 gen23_{2} + 0.3333333333333333 gen23_{3} + 0.3333333333333333 gen23_{4} $"
],
"text/plain": [
"0.3333333333333333 Pgen21 + 0.3333333333333333 gen21[1] + 0.3333333333333333 gen21[2] + 0.3333333333333333 gen21[3] + 0.3333333333333333 gen21[4] + 0.3333333333333333 Pgen22 + 0.3333333333333333 gen22[1] + 0.3333333333333333 gen22[2] + 0.3333333333333333 gen22[3] + 0.3333333333333333 gen22[4] + 0.3333333333333333 Pgen23 + 0.3333333333333333 gen23[1] + 0.3333333333333333 gen23[2] + 0.3333333333333333 gen23[3] + 0.3333333333333333 gen23[4]"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cost2 = [data2[s][\"cost\"] for s=1:ns]\n",
"cost2_exp = cost2' * proba_scenar"
]
},
{
"cell_type": "markdown",
"id": "e617d88b-c5a7-4c38-aef0-ebd71b92dbfc",
"metadata": {},
"source": [
"Coupling between the two stages"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "b0383822-24dc-4b8d-a79e-7ca9f26ed5a8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:\n",
" Pgen1 - Pgen21 ≤ 0\n",
" Pgen1 - Pgen22 ≤ 0\n",
" Pgen1 - Pgen23 ≤ 0"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@constraint(m, [s=1:ns], data1[\"Pgen\"] <= data2[s][\"Pgen\"])"
]
},
{
"cell_type": "markdown",
"id": "7102487b-aa50-4a53-bfda-4507d3c46400",
"metadata": {},
"source": [
"Set optimization objective, summing both stages"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "f26c6823-06ae-4722-9d65-bf3fb7811a80",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"A JuMP Model\n",
"Minimization problem with:\n",
"Variables: 20\n",
"Objective function type: AffExpr\n",
"`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 16 constraints\n",
"`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 19 constraints\n",
"Model mode: AUTOMATIC\n",
"CachingOptimizer state: EMPTY_OPTIMIZER\n",
"Solver name: HiGHS"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@objective(m, Min, data1[\"cost\"] + cost2_exp)\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "e197a786-6ff3-4931-b984-9bfe8b8df1f2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pgen1=1.5\n",
"Pgen2=[1.5, 1.5, 3.0]\n"
]
}
],
"source": [
"optimize!(m)\n",
"Pgen1_val = value(data1[\"Pgen\"])\n",
"println(\"Pgen1=$(Pgen1_val)\")\n",
"Pgen2_val = [value(data2[s][\"Pgen\"]) for s=1:ns]\n",
"println(\"Pgen2=$(Pgen2_val)\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.10.4",
"language": "julia",
"name": "julia-1.10"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
@pierre-haessig
Copy link
Author

Associated post on Julia forum: https://discourse.julialang.org/t/modular-implementation-of-optimization-problems-with-jump/116261

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment