{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Words of warning"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"SDDP is a powerful solution technique for multistage stochastic programming.\n",
"However, there are a number of subtle things to be aware of before creating\n",
"your own models."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Relatively complete recourse"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Models built in SDDP.jl need a property called _relatively complete recourse_."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"One definition of relatively complete recourse is that _all_ feasible decisions\n",
"(not necessarily optimal) in a subproblem lead to feasible decisions in future\n",
"subproblems."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"For example, in the following problem, one feasible first stage decision is\n",
"`x.out = 0`. But this causes an infeasibility in the second stage which requires\n",
"`x.in >= 1`. This will throw an error about infeasibility if you try to solve."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using SDDP, HiGHS\n",
"\n",
"model = SDDP.LinearPolicyGraph(\n",
" stages = 2,\n",
" lower_bound = 0,\n",
" optimizer = HiGHS.Optimizer,\n",
") do sp, t\n",
" @variable(sp, x >= 0, SDDP.State, initial_value = 1)\n",
" if t == 2\n",
" @constraint(sp, x.in >= 1)\n",
" end\n",
" @stageobjective(sp, x.out)\n",
"end\n",
"\n",
"try #hide\n",
" SDDP.train(model, iteration_limit = 1, print_level = 0)\n",
"catch err #hide\n",
" showerror(stderr, err) #hide\n",
"end #hide"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"!!! warning\n",
" The actual constraints causing the infeasibilities can be deceptive! A good\n",
" strategy to debug is to comment out all constraints. Then, one-by-one,\n",
" un-comment the constraints and try resolving the model to check if it finds a\n",
" feasible solution."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Numerical stability"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"If you aren't aware, SDDP builds an outer-approximation to a convex function\n",
"using cutting planes. This results in a formulation that is particularly hard\n",
"for solvers like HiGHS, Gurobi, and CPLEX to deal with. As a result, you may\n",
"run into weird behavior. This behavior could include:\n",
"\n",
" - Iterations suddenly taking a long time (the solver stalled)\n",
" - Subproblems turning infeasible or unbounded after many iterations\n",
" - Solvers returning \"Numerical Error\" statuses"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Problem scaling"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"In almost all cases, the cause of this is poor problem scaling. For our\n",
"purpose, poor problem scaling means having variables with very large numbers\n",
"and variables with very small numbers in the same model."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"!!! tip\n",
" Gurobi has an excellent [set of articles](http://www.gurobi.com/documentation/8.1/refman/numerics_gurobi_guidelines.html)\n",
" on numerical issues and how to avoid them."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Consider, for example, the hydro-thermal scheduling problem we have been\n",
"discussing in previous tutorials."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"If we define the volume of the reservoir in terms of m³, then a lake might\n",
"have a capacity of 10^10 m³: `@variable(subproblem, 0 <= volume <= 10^10)`.\n",
"Moreover, the cost per cubic meter might be around \\$0.05/m³. To calculate\n",
"the value of water in our reservoir, we need to multiple a variable on the\n",
"order of 10^10, by one on the order of 10⁻²! That is twelve orders of\n",
"magnitude!"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"To improve the performance of the SDDP algorithm (and reduce the chance of\n",
"weird behavior), try to re-scale the units of the problem in order to reduce\n",
"the largest difference in magnitude. For example, if we talk in terms of\n",
"million m³, then we have a capacity of 10⁴ million m³, and a price of\n",
"\\$50,000 per million m³. Now things are only one order of magnitude apart."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Numerical stability report"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"To aid in the diagnose of numerical issues, you can call\n",
"`SDDP.numerical_stability_report`. By default, this aggregates all of\n",
"the nodes into a single report. You can produce a stability report for each\n",
"node by passing `by_node=true`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using SDDP\n",
"\n",
"model = SDDP.LinearPolicyGraph(stages = 2, lower_bound = -1e10) do subproblem, t\n",
" @variable(subproblem, x >= -1e7, SDDP.State, initial_value = 1e-5)\n",
" @constraint(subproblem, 1e9 * x.out >= 1e-6 * x.in + 1e-8)\n",
" @stageobjective(subproblem, 1e9 * x.out)\n",
"end\n",
"\n",
"SDDP.numerical_stability_report(model)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The report analyses the magnitude (in absolute terms) of the coefficients in\n",
"the constraint matrix, the objective function, any variable bounds, and in the\n",
"RHS of the constraints. A warning will be thrown in `SDDP.jl` detects very\n",
"large or small values. As discussed in Problem scaling, this is an\n",
"indication that you should reformulate your model."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"By default, a numerical stability check is run when you call\n",
"`SDDP.train`, although it can be turned off by passing\n",
"`run_numerical_stability_report = false`."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Solver-specific options"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"If you have a particularly troublesome model, you should investigate setting\n",
"solver-specific options to improve the numerical stability of each solver. For\n",
"example, Gurobi has a [`NumericFocus`\n",
"option](http://www.gurobi.com/documentation/8.1/refman/numericfocus.html#parameter:NumericFocus)."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Choosing an initial bound"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"One of the important requirements when building a SDDP model is to choose an\n",
"appropriate bound on the objective (lower if minimizing, upper if maximizing).\n",
"However, it can be hard to choose a bound if you don't know the solution!\n",
"(Which is very likely.)"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"The bound should not be as large as possible (since this will help with\n",
"convergence and the numerical issues discussed above), but if chosen too\n",
"small, it may cut off the feasible region and lead to a sub-optimal solution."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Consider the following simple model, where we first set `lower_bound` to `0.0`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using SDDP, HiGHS\n",
"\n",
"model = SDDP.LinearPolicyGraph(\n",
" stages = 3,\n",
" sense = :Min,\n",
" lower_bound = 0.0,\n",
" optimizer = HiGHS.Optimizer,\n",
") do subproblem, t\n",
" @variable(subproblem, x >= 0, SDDP.State, initial_value = 2)\n",
" @variable(subproblem, u >= 0)\n",
" @variable(subproblem, v >= 0)\n",
" @constraint(subproblem, x.out == x.in - u)\n",
" @constraint(subproblem, u + v == 1.5)\n",
" @stageobjective(subproblem, t * v)\n",
"end\n",
"\n",
"SDDP.train(model, iteration_limit = 5, run_numerical_stability_report = false)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"Now consider the case when we set the `lower_bound` to `10.0`:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using SDDP, HiGHS\n",
"\n",
"model = SDDP.LinearPolicyGraph(\n",
" stages = 3,\n",
" sense = :Min,\n",
" lower_bound = 10.0,\n",
" optimizer = HiGHS.Optimizer,\n",
") do subproblem, t\n",
" @variable(subproblem, x >= 0, SDDP.State, initial_value = 2)\n",
" @variable(subproblem, u >= 0)\n",
" @variable(subproblem, v >= 0)\n",
" @constraint(subproblem, x.out == x.in - u)\n",
" @constraint(subproblem, u + v == 1.5)\n",
" @stageobjective(subproblem, t * v)\n",
"end\n",
"\n",
"SDDP.train(model, iteration_limit = 5, run_numerical_stability_report = false)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"How do we tell which is more appropriate? There are a few clues that you\n",
"should look out for.\n",
"\n",
"- The bound converges to a value above (if minimizing) the simulated cost of\n",
" the policy. In this case, the problem is deterministic, so it is easy to\n",
" tell. But you can also check by performing a Monte Carlo simulation like we\n",
" did in An introduction to SDDP.jl.\n",
"\n",
"- The bound converges to different values when we change the bound. This is\n",
" another clear give-away. The bound provided by the user is only used in the\n",
" initial iterations. __It should not change the value of the converged\n",
" policy.__ Thus, if you don't know an appropriate value for the bound, choose\n",
" an initial value, and then increase (or decrease) the value of the bound to\n",
" confirm that the value of the policy doesn't change.\n",
"\n",
"- The bound converges to a value _close_ to the bound provided by the user.\n",
" This varies between models, but notice that `11.0` is quite close to `10.0`\n",
" compared with `3.5` and `0.0`."
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.4"
},
"kernelspec": {
"name": "julia-1.10",
"display_name": "Julia 1.10.4",
"language": "julia"
}
},
"nbformat": 4
}