diff --git a/Project.toml b/Project.toml index 9c3832e..6aa6aeb 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,8 @@ PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" [compat] DataStructures = "0.18" @@ -26,5 +28,7 @@ GZip = "0.5" JSON = "0.21" JuMP = "1" MathOptInterface = "1" +MPI = "0.20.8" PackageCompiler = "1" julia = "1" +TimerOutputs = "0.5.23" \ No newline at end of file diff --git a/docs/src/format.md b/docs/src/format.md index 267da58..24b8de3 100644 --- a/docs/src/format.md +++ b/docs/src/format.md @@ -4,7 +4,7 @@ Data Format Input Data Format ----------------- -Instances are specified by JSON files containing the following main sections: +To create the JuMP model for a two-stage security constrained unit commitment (SUC) problem, UC.jl creates an instance that is composed of one or multiple scenarios. For each scenario, UC.jl expects a JSON file that describes the parameters of the problem, as well as the weight and the name of each respective scenario. Note that the input data format required for modeling the deterministic security-constrained unit commitment (SCUC) problem also follows the format mapped out below, as SCUC is treated by the package as the special case of the SUC problem with one scenario. The requisite input data format contains the following fields: * [Parameters](#Parameters) * [Buses](#Buses) @@ -26,6 +26,8 @@ This section describes system-wide parameters, such as power balance penalty, an | `Time horizon (h)` | Length of the planning horizon (in hours). | Required | N | `Time step (min)` | Length of each time step (in minutes). Must be a divisor of 60 (e.g. 60, 30, 20, 15, etc). | `60` | N | `Power balance penalty ($/MW)` | Penalty for system-wide shortage or surplus in production (in $/MW). This is charged per time step. For example, if there is a shortage of 1 MW for three time steps, three times this amount will be charged. | `1000.0` | Y +| `Scenario name` | Name of the scenario. | `s1` | N +| `Scenario weight` | Weight of the scenario. The scenario weight can be any positive real number, that is, it does not have to be between zero and one. The package normalizes the weights to ensure that the probability of all scenarios sum up to one. | one divided by the total number of scenarios | N #### Example @@ -34,7 +36,9 @@ This section describes system-wide parameters, such as power balance penalty, an "Parameters": { "Version": "0.3", "Time horizon (h)": 4, - "Power balance penalty ($/MW)": 1000.0 + "Power balance penalty ($/MW)": 1000.0, + "Scenario name": "s1", + "Scenario weight": 0.5 } } ``` diff --git a/docs/src/index.md b/docs/src/index.md index 9a8a013..752fb64 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,6 +1,6 @@ # UnitCommitment.jl -**UnitCommitment.jl** (UC.jl) is a Julia/JuMP optimization package for the Security-Constrained Unit Commitment Problem (SCUC), a fundamental optimization problem in power systems used, for example, to clear the day-ahead electricity markets. The package provides benchmark instances for the problem and Julia/JuMP implementations of state-of-the-art mixed-integer programming formulations. +**UnitCommitment.jl** (UC.jl) is a Julia/JuMP optimization package for modeling and solving the two-stage Stochastic Unit Commitment (SUC) problem, which permits the representation of uncertainty in a broad range of problem parameters. As a special case of the SUC problem, UC.jl allows for modeling and solving the Security-Constrained Unit Commitment Problem (SCUC), a fundamental optimization problem in power systems used, for example, to clear the day-ahead electricity markets. The package provides benchmark instances for the problem and Julia/JuMP implementations of state-of-the-art mixed-integer programming formulations. ## Package Components diff --git a/docs/src/model.md b/docs/src/model.md index d944b07..f57648d 100644 --- a/docs/src/model.md +++ b/docs/src/model.md @@ -5,9 +5,16 @@ In this page, we describe the JuMP optimization model produced by the function ` Decision variables ------------------ +In this section, we describe the decision variables of the JuMP model produced by the function `UnitCommitment.build_model`. + +Recall that the UC.jl package allows for modeling and solving the SUC problem, wherein a broad range problem parameters can vary across user-defined scenarios. To facilitate the modeling of the SUC problem, some of the decision variables of the JuMP model are modeled as first-stage decisions, which are here-and-now decisions taken before the uncertainty is realized and thus are the same across all scenarios. As laid out below, the first-stage decisions relate primarily to the binary decisions of the generators. + +In contrast, most of the decision variables of the JuMP model are second-stage decisions that can attain a different value in each scenario. Accordingly, the second-stage decision variables described below are indexed by the term `sn`. The term `sn` is preserved even if the package is used to model and solve the deterministic security-constrained unit commitment (SCUC) problem, which is currently treated by the package as the special case of the SUC model with only one scenario, having the default scenario index `"s1"`. ### Generators +In this section, we describe the decision variables associated with the generators, which include both thermal units (e.g., natural gas-fired power plant) and profiled units (e.g., wind turbine). + #### Thermal Units Name | Symbol | Description | Unit @@ -15,11 +22,15 @@ Name | Symbol | Description | Unit `is_on[g,t]` | $u_{g}(t)$ | True if generator `g` is on at time `t`. | Binary `switch_on[g,t]` | $v_{g}(t)$ | True is generator `g` switches on at time `t`. | Binary `switch_off[g,t]` | $w_{g}(t)$ | True if generator `g` switches off at time `t`. | Binary -`prod_above[g,t]` |$p'_{g}(t)$ | Amount of power produced by generator `g` above its minimum power output at time `t`. For example, if the minimum power of generator `g` is 100 MW and `g` is producing 115 MW of power at time `t`, then `prod_above[g,t]` equals `15.0`. | MW -`segprod[g,t,k]` | $p^k_g(t)$ | Amount of power from piecewise linear segment `k` produced by generator `g` at time `t`. For example, if cost curve for generator `g` is defined by the points `(100, 1400)`, `(110, 1600)`, `(130, 2200)` and `(135, 2400)`, and if the generator is producing 115 MW of power at time `t`, then `segprod[g,t,:]` equals `[10.0, 5.0, 0.0]`.| MW -`reserve[r,g,t]` | $r_g(t)$ | Amount of reserve `r` provided by unit `g` at time `t`. | MW `startup[g,t,s]` | $\delta^s_g(t)$ | True if generator `g` switches on at time `t` incurring start-up costs from start-up category `s`. | Binary +`prod_above[sn,g,t]` |$p'_{g}(t)$ | Amount of power produced by generator `g` above its minimum power output at time `t` in scenario `sn`. For example, if the minimum power of generator `g` is 100 MW and `g` is producing 115 MW of power at time `t` in scenario `sn`, then `prod_above[sn,g,t]` equals `15.0`. | MW +`segprod[sn,g,t,k]` | $p^k_g(t)$ | Amount of power from piecewise linear segment `k` produced by generator `g` at time `t` in scenario `sn`. For example, if cost curve for generator `g` is defined by the points `(100, 1400)`, `(110, 1600)`, `(130, 2200)` and `(135, 2400)`, and if the generator is producing 115 MW of power at time `t` in scenario `sn`, then `segprod[sn,g,t,:]` equals `[10.0, 5.0, 0.0]`.| MW +`reserve[sn,r,g,t]` | $r_g(t)$ | Amount of reserve `r` provided by unit `g` at time `t` in scenario `sn`. | MW + +!!! warning + The first-stage decision variables of the JuMP model are `is_on[g,t]`, `switch_on[g,t]`, `switch_off[g,t]`, and `startup[g,t,s]`. As such, the dictionaries corresponding to these variables do not include the scenario index in their keys. In contrast, all other variables of the created JuMP model are allowed to obtain a different value in each scenario and are thus modeled as second-stage decision variables. Accordingly, the dictionaries of all second-stage decision variables have the scenario index in their keys. This is true even if the model is created to solve the deterministic SCUC, in which case the default scenario index `s1` is included in the dictionary key. + #### Profiled Units @@ -32,26 +43,26 @@ Name | Symbol | Description | Unit Name | Symbol | Description | Unit :-----|:------:|:-------------|:------: -`net_injection[b,t]` | $n_b(t)$ | Net injection at bus `b` at time `t`. | MW -`curtail[b,t]` | $s^+_b(t)$ | Amount of load curtailed at bus `b` at time `t` | MW +`net_injection[sn,b,t]` | $n_b(t)$ | Net injection at bus `b` at time `t` in scenario `sn`. | MW +`curtail[sn,b,t]` | $s^+_b(t)$ | Amount of load curtailed at bus `b` at time `t` in scenario `sn`. | MW ### Price-sensitive loads Name | Symbol | Description | Unit :-----|:------:|:-------------|:------: -`loads[s,t]` | $d_{s}(t)$ | Amount of power served to price-sensitive load `s` at time `t`. | MW +`loads[sn,s,t]` | $d_{s}(t)$ | Amount of power served to price-sensitive load `s` at time `t` in scenario `sn`. | MW ### Transmission lines Name | Symbol | Description | Unit :-----|:------:|:-------------|:------: -`flow[l,t]` | $f_l(t)$ | Power flow on line `l` at time `t`. | MW -`overflow[l,t]` | $f^+_l(t)$ | Amount of flow above the limit for line `l` at time `t`. | MW +`flow[sn,l,t]` | $f_l(t)$ | Power flow on line `l` at time `t` in scenario `sn`. | MW +`overflow[sn,l,t]` | $f^+_l(t)$ | Amount of flow above the limit for line `l` at time `t` in scenario `sn`. | MW !!! warning - Since transmission and N-1 security constraints are enforced in a lazy way, most of the `flow[l,t]` variables are never added to the model. Accessing `model[:flow][l,t]` without first checking that the variable exists will likely generate an error. + Since transmission and N-1 security constraints are enforced in a lazy way, most of the `flow[l,t]` variables are never added to the model. Accessing `model[:flow][sn,l,t]` without first checking that the variable exists will likely generate an error. Objective function ------------------ @@ -106,7 +117,12 @@ end ### Fixing variables, modifying objective function and adding constraints -Since we now have a direct reference to the JuMP decision variables, it is possible to fix variables, change the coefficients in the objective function, or even add new constraints to the model before solving it. The script below shows how can this be accomplished. For more information on modifying an existing model, [see the JuMP documentation](https://jump.dev/JuMP.jl/stable/manual/variables/). +Since we now have a direct reference to the JuMP decision variables, it is possible to fix variables, change the coefficients in the objective function, or even add new constraints to the model before solving it. +!!! warning + + It is important to take into account the stage of the decision variables in modifying the optimization model. In changing a deterministic SCUC model, modifying the second-stage decision variables requires adding the term `s1`, which is the default scenario name assigned to the second-stage decision variables in the SCUC model. For an SUC model, the package permits the modification of the second-stage decision variables individually for each scenario. + +The script below shows how the JuMP model can be modified after it is created. For more information on modifying an existing model, [see the JuMP documentation](https://jump.dev/JuMP.jl/stable/manual/variables/). ```julia using Cbc @@ -122,13 +138,31 @@ model = UnitCommitment.build_model( optimizer=Cbc.Optimizer, ) -# Fix a decision variable to 1.0 +# Fix the commitment status of the generator "g1" in time period 1 to 1.0 JuMP.fix( model[:is_on]["g1",1], 1.0, force=true, ) +# Fix the production level of the generator "g1" above its minimum level in time period 1 and +# in scenario "s1" to 20.0 MW. Observe that the three-tuple dictionary key involves the scenario +# index "s1", as production above minimum is a second-stage decision variable. + +JuMP.fix( + model[:prod_above]["s1", "g1", 1], + 20.0, + force=true, +) + +# Enforce the curtailment of 20.0 MW of load at bus "b2" in time period 4 in scenario "s10". +JuMP.fix( + curtail["s10", "b2", 4] = + 20.0, + force=true, +) + + # Change the objective function JuMP.set_objective_coefficient( model, @@ -178,10 +212,10 @@ for t in 1:T # In this example, we assume a cost of $5/MW. set_objective_coefficient(model, x[t], 5.0) - # Attach the new component to bus b1, by modifying the + # Attach the new component to bus b1 in scenario s1, by modifying the # constraint `eq_net_injection`. set_normalized_coefficient( - model[:eq_net_injection]["b1", t], + model[:eq_net_injection]["s1", "b1", t], x[t], 1.0, ) diff --git a/docs/src/usage.md b/docs/src/usage.md index 26d1779..cc0281a 100644 --- a/docs/src/usage.md +++ b/docs/src/usage.md @@ -25,7 +25,7 @@ Typical Usage ### Solving user-provided instances -The first step to use UC.jl is to construct a JSON file describing your unit commitment instance. See [Data Format](format.md) for a complete description of the data format UC.jl expects. The next steps, as shown below, are to: (1) read the instance from file; (2) construct the optimization model; (3) run the optimization; and (4) extract the optimal solution. +The first step to use UC.jl is to construct a JSON file that describes each scenario of your unit commitment instance. See [Data Format](format.md) for a complete description of the data format UC.jl expects. The next steps, as shown below, are to: (1) construct the instance using scenario files; (2) build the optimization model; (3) run the optimization; and (4) extract the optimal solution. ```julia using Cbc @@ -33,7 +33,7 @@ using JSON using UnitCommitment # 1. Read instance -instance = UnitCommitment.read("/path/to/input.json") +instance = UnitCommitment.read(["/path/to/s1.json", "/path/to/s2.json"]) # 2. Construct optimization model model = UnitCommitment.build_model( @@ -49,6 +49,15 @@ solution = UnitCommitment.solution(model) UnitCommitment.write("/path/to/output.json", solution) ``` +The above lines of code can also be used for solving the deterministic Security-Constrained Unit Commitment (SCUC) problem, which will create an instance based on a single scenario. The unit commitment instance for SCUC can alternatively be constructed as + +```julia +# 1. Read instance +instance = UnitCommitment.read("/path/to/input.json") +``` + + + ### Solving benchmark instances UnitCommitment.jl contains a large number of benchmark instances collected from the literature and converted into a common data format. To solve one of these instances individually, instead of constructing your own, the function `read_benchmark` can be used, as shown below. See [Instances](instances.md) for the complete list of available instances. @@ -136,6 +145,80 @@ solution = JSON.parsefile("solution.json") # Validate solution and print validation errors UnitCommitment.validate(instance, solution) ``` +## Modeling and Solving the SUC Problem + +To model the SUC problem, UC.jl supports reading scenario files at a specified directory using the `Glob` package. For instance, in order to construct a SUC instance using all JSON files at a given directory, where each JSON file describes one scenario, the following line of code can be used: + +```julia +instance = UnitCommitment.read(glob("*.json", "/path/to/scenarios/")) +``` + +Alternatively, the specific vector of scenario files can also be passed as follows: +```julia +instance = UnitCommitment.read(["/path/to/s1.json", "/path/to/s2.json"])) +``` + +## Solving the SUC Problem + +We next lay out the alternative methods supported by UC.jl for solving the SUC problem. + +## Solving the Extensive Form of the SUC Problem + +By default, UC.jl solves the extensive form of the SUC problem. + +```julia +UnitCommitment.optimize!(model) +solution = UnitCommitment.solution(model) +UnitCommitment.write("/path/to/output.json", solution) +``` + +Note that the created `solution` dictionary will include both the optimal first-stage decisions, as well as the optimal second-stage decisions under all scenarios. + +## Solving the SUC Problem Using Progressive Hedging + +Importantly, UC.jl further provides the option of solving the SUC problem using the progressive hedging (PH) algorithm, which is an algorithm closely related to the alternating direction method of multipliers (ADMM). To that end, the package supports solving the PH subproblem associated with each scenario in parallel in a separate Julia process, where the communication among the Julia processes is provided using the Message Passing Interface or MPI. + +The solve the SUC problem using Progressive Hedging, you may run the following line of code, which will create `NUM_OF_PROCS` processes where each process executes the `ph.jl` file. + +```julia +using MPI +const FILENAME = "ph.jl" +const NUM_OF_PROCS = 5 + +mpiexec(exe -> run(`$exe -n $NUM_OF_PROCS $(Base.julia_cmd()) $FILENAME`)) +``` + +#### **`ph.jl`** +```julia +using MPI: MPI_Info +using Gurobi, MPI, UnitCommitment, Glob +import UnitCommitment: ProgressiveHedging + +MPI.Init() +ph = ProgressiveHedging.Method() +instance = UnitCommitment.read( + glob("*.json", "/path/to/scenarios/"), + ph +) +model = UnitCommitment.build_model( + instance = instance, + optimizer = Gurobi.Optimizer, + ) +UnitCommitment.optimize!(model, ph) +solution = UnitCommitment.solution(model, ph) +MPI.Finalize() +``` + +Observe that the `read`, `build_model`, and `solution` methods take the `ph` object as an argument, which is of type `Progressive Hedging`. + +The subproblem solved within each Julia process deduces the number of scenarios it needs to model using the total number of scenarios and the total number of processes. For instance, if `glob("*.json", "/path/to/scenarios/")` returns a vector of 15 scenario file paths and `NUM_OF_PROCS = 5`, then each subproblem will model and solve 3 scenarios. If the total number of scenarios is not divisible by `NUM_OF_PROCS`, then the read method will throw an error. + +The `solution(model, ph)` method gathers the optimal solution of all processes and returns a dictionary that contains all optimal first-stage decisions as well as all optimal second-stage decisions evaluated for each scenario. + +!!! warning + + Currently, PH can handle only equiprobable scenarios. Further, `solution(model, ph)` can only handle cases where only one scenario is modeled in each process. + ## Computing Locational Marginal Prices diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 3191395..40016f7 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -19,6 +19,7 @@ include("model/formulations/KnuOstWat2018/structs.jl") include("model/formulations/MorLatRam2013/structs.jl") include("model/formulations/PanGua2016/structs.jl") include("solution/methods/XavQiuWanThi2019/structs.jl") +include("solution/methods/ProgressiveHedging/structs.jl") include("model/formulations/WanHob2016/structs.jl") include("import/egret.jl") @@ -49,6 +50,9 @@ include("solution/methods/XavQiuWanThi2019/enforce.jl") include("solution/methods/XavQiuWanThi2019/filter.jl") include("solution/methods/XavQiuWanThi2019/find.jl") include("solution/methods/XavQiuWanThi2019/optimize.jl") +include("solution/methods/ProgressiveHedging/optimize.jl") +include("solution/methods/ProgressiveHedging/read.jl") +include("solution/methods/ProgressiveHedging/solution.jl") include("solution/optimize.jl") include("solution/solution.jl") include("solution/warmstart.jl") diff --git a/src/solution/methods/ProgressiveHedging/optimize.jl b/src/solution/methods/ProgressiveHedging/optimize.jl new file mode 100644 index 0000000..c6831e5 --- /dev/null +++ b/src/solution/methods/ProgressiveHedging/optimize.jl @@ -0,0 +1,265 @@ +# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. +using MPI, Printf +using TimerOutputs +import JuMP +const to = TimerOutput() + +function optimize!( + model::JuMP.Model, + method::ProgressiveHedging.Method, +)::ProgressiveHedging.FinalResult + mpi = ProgressiveHedging.MpiInfo(MPI.COMM_WORLD) + iterations = Array{ProgressiveHedging.IterationInfo,1}(undef, 0) + if method.consensus_vars === nothing + method.consensus_vars = + [var for var in all_variables(model) if is_binary(var)] + end + nvars = length(method.consensus_vars) + if method.weights === nothing + method.weights = [1.0 for _ in 1:nvars] + end + if method.initial_global_consensus_vals === nothing + method.initial_global_consensus_vals = [0.0 for _ in 1:nvars] + end + + ph_sp_params = ProgressiveHedging.SpParams( + ρ = method.ρ, + λ = [method.λ_default for _ in 1:nvars], + global_consensus_vals = method.initial_global_consensus_vals, + ) + ph_subproblem = ProgressiveHedging.SubProblem( + model, + model[:obj], + method.consensus_vars, + method.weights, + ) + set_optimizer_attribute(model, "Threads", method.num_of_threads) + while true + it_time = @elapsed begin + solution = solve_subproblem(ph_subproblem, ph_sp_params) + MPI.Barrier(mpi.comm) + global_obj = compute_global_objective(mpi, solution) + global_consensus_vals = compute_global_consensus(mpi, solution) + update_λ_and_residuals!( + solution, + ph_sp_params, + global_consensus_vals, + ) + global_infeas = compute_global_infeasibility(solution, mpi) + global_residual = compute_global_residual(mpi, solution) + if has_numerical_issues(global_consensus_vals) + break + end + end + total_elapsed_time = compute_total_elapsed_time(it_time, iterations) + it = ProgressiveHedging.IterationInfo( + it_num = length(iterations) + 1, + sp_consensus_vals = solution.consensus_vals, + global_consensus_vals = global_consensus_vals, + sp_obj = solution.obj, + global_obj = global_obj, + it_time = it_time, + total_elapsed_time = total_elapsed_time, + global_residual = global_residual, + global_infeas = global_infeas, + ) + iterations = [iterations; it] + print_progress(mpi, it, method.print_interval) + if should_stop(mpi, iterations, method.termination_criteria) + break + end + end + + return ProgressiveHedging.FinalResult( + last(iterations).global_obj, + last(iterations).sp_consensus_vals, + last(iterations).global_infeas, + last(iterations).it_num, + last(iterations).total_elapsed_time, + ) +end + +function compute_total_elapsed_time( + it_time::Float64, + iterations::Array{ProgressiveHedging.IterationInfo,1}, +)::Float64 + length(iterations) > 0 ? + current_total_time = last(iterations).total_elapsed_time : + current_total_time = 0 + return current_total_time + it_time +end + +function compute_global_objective( + mpi::ProgressiveHedging.MpiInfo, + s::ProgressiveHedging.SpSolution, +)::Float64 + global_obj = MPI.Allreduce(s.obj, MPI.SUM, mpi.comm) + global_obj /= mpi.nprocs + return global_obj +end + +function compute_global_consensus( + mpi::ProgressiveHedging.MpiInfo, + s::ProgressiveHedging.SpSolution, +)::Array{Float64,1} + sp_consensus_vals = s.consensus_vals + global_consensus_vals = MPI.Allreduce(sp_consensus_vals, MPI.SUM, mpi.comm) + global_consensus_vals = global_consensus_vals / mpi.nprocs + return global_consensus_vals +end + +function compute_global_residual( + mpi::ProgressiveHedging.MpiInfo, + s::ProgressiveHedging.SpSolution, +)::Float64 + n_vars = length(s.consensus_vals) + local_residual_sum = abs.(s.residuals) + global_residual_sum = MPI.Allreduce(local_residual_sum, MPI.SUM, mpi.comm) + return sum(global_residual_sum) / n_vars +end + +function compute_global_infeasibility( + solution::ProgressiveHedging.SpSolution, + mpi::ProgressiveHedging.MpiInfo, +)::Float64 + local_infeasibility = norm(solution.residuals) + global_infeas = MPI.Allreduce(local_infeasibility, MPI.SUM, mpi.comm) + return global_infeas +end + +function solve_subproblem( + sp::ProgressiveHedging.SubProblem, + ph_sp_params::ProgressiveHedging.SpParams, +)::ProgressiveHedging.SpSolution + G = length(sp.consensus_vars) + if norm(ph_sp_params.λ) < 1e-3 + @objective(sp.mip, Min, sp.obj) + else + @objective( + sp.mip, + Min, + sp.obj + + sum( + sp.weights[g] * + ph_sp_params.λ[g] * + (sp.consensus_vars[g] - ph_sp_params.global_consensus_vals[g]) + for g in 1:G + ) + + (ph_sp_params.ρ / 2) * sum( + sp.weights[g] * + ( + sp.consensus_vars[g] - + ph_sp_params.global_consensus_vals[g] + )^2 for g in 1:G + ) + ) + end + optimize!(sp.mip, XavQiuWanThi2019.Method()) + obj = objective_value(sp.mip) + sp_consensus_vals = value.(sp.consensus_vars) + return ProgressiveHedging.SpSolution( + obj = obj, + consensus_vals = sp_consensus_vals, + residuals = zeros(G), + ) +end + +function update_λ_and_residuals!( + solution::ProgressiveHedging.SpSolution, + ph_sp_params::ProgressiveHedging.SpParams, + global_consensus_vals::Array{Float64,1}, +)::Nothing + n_vars = length(solution.consensus_vals) + ph_sp_params.global_consensus_vals = global_consensus_vals + for n in 1:n_vars + solution.residuals[n] = + solution.consensus_vals[n] - ph_sp_params.global_consensus_vals[n] + ph_sp_params.λ[n] += ph_sp_params.ρ * solution.residuals[n] + end +end + +function print_header(mpi::ProgressiveHedging.MpiInfo)::Nothing + if !mpi.root + return + end + @info "Solving via Progressive Hedging:" + @info @sprintf( + "%8s %20s %20s %14s %8s %8s", + "iter", + "obj", + "infeas", + "consensus", + "time-it", + "time" + ) +end + +function print_progress( + mpi::ProgressiveHedging.MpiInfo, + iteration::ProgressiveHedging.IterationInfo, + print_interval, +)::Nothing + if !mpi.root + return + end + if iteration.it_num % print_interval != 0 + return + end + @info @sprintf( + "Current iteration %8d %20.6e %20.6e %12.2f %% %8.2f %8.2f", + iteration.it_num, + iteration.global_obj, + iteration.global_infeas, + iteration.global_residual * 100, + iteration.it_time, + iteration.total_elapsed_time + ) +end + +function has_numerical_issues(target::Array{Float64,1})::Bool + if target == NaN + @warn "Numerical issues detected. Stopping." + return true + end + return false +end + +function should_stop( + mpi::ProgressiveHedging.MpiInfo, + iterations::Array{ProgressiveHedging.IterationInfo,1}, + criteria::ProgressiveHedging.TerminationCriteria, +)::Bool + if length(iterations) >= criteria.max_iterations + if mpi.root + @info "Iteration limit reached. Stopping." + end + return true + end + + if length(iterations) < criteria.min_iterations + return false + end + + if last(iterations).total_elapsed_time > criteria.max_time + if mpi.root + @info "Time limit reached. Stopping." + end + return true + end + + curr_it = last(iterations) + prev_it = iterations[length(iterations)-1] + + if curr_it.global_infeas < criteria.min_feasibility + obj_change = abs(prev_it.global_obj - curr_it.global_obj) + if obj_change < criteria.min_improvement + if mpi.root + @info "Feasibility limit reached. Stopping." + end + return true + end + end + return false +end diff --git a/src/solution/methods/ProgressiveHedging/read.jl b/src/solution/methods/ProgressiveHedging/read.jl new file mode 100644 index 0000000..476b8b7 --- /dev/null +++ b/src/solution/methods/ProgressiveHedging/read.jl @@ -0,0 +1,18 @@ +# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +function read( + paths::Vector{String}, + method::ProgressiveHedging.Method, +)::UnitCommitmentInstance + comm = MPI.COMM_WORLD + mpi = ProgressiveHedging.MpiInfo(comm) + (length(paths) % mpi.nprocs == 0) || error( + "Number of processes $(mpi.nprocs) is not a divisor of $(length(paths))", + ) + bundled_scenarios = length(paths) ÷ mpi.nprocs + sc_num_start = (mpi.rank - 1) * bundled_scenarios + 1 + sc_num_end = mpi.rank * bundled_scenarios + return read(paths[sc_num_start:sc_num_end]) +end diff --git a/src/solution/methods/ProgressiveHedging/solution.jl b/src/solution/methods/ProgressiveHedging/solution.jl new file mode 100644 index 0000000..b06ad89 --- /dev/null +++ b/src/solution/methods/ProgressiveHedging/solution.jl @@ -0,0 +1,86 @@ +# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using MPI, DataStructures +const FIRST_STAGE_VARS = ["Is on", "Switch on", "Switch off"] + +function solution( + model::JuMP.Model, + method::ProgressiveHedging.Method, +)::OrderedDict + comm = MPI.COMM_WORLD + mpi = ProgressiveHedging.MpiInfo(comm) + sp_solution = UnitCommitment.solution(model) + gather_solution = OrderedDict() + for (solution_key, dict) in sp_solution + if solution_key !== "Spinning reserve (MW)" && + solution_key ∉ FIRST_STAGE_VARS + push!(gather_solution, solution_key => OrderedDict()) + for (gen_bus_key, values) in dict + global T = length(values) + receive_values = + MPI.UBuffer(Vector{Float64}(undef, T * mpi.nprocs), T) + MPI.Gather!(float.(values), receive_values, comm) + if mpi.root + push!( + gather_solution[solution_key], + gen_bus_key => receive_values.data, + ) + end + end + end + end + push!(gather_solution, "Spinning reserve (MW)" => OrderedDict()) + for (reserve_type, dict) in sp_solution["Spinning reserve (MW)"] + push!( + gather_solution["Spinning reserve (MW)"], + reserve_type => OrderedDict(), + ) + for (gen_key, values) in dict + receive_values = + MPI.UBuffer(Vector{Float64}(undef, T * mpi.nprocs), T) + MPI.Gather!(float.(values), receive_values, comm) + if mpi.root + push!( + gather_solution["Spinning reserve (MW)"][reserve_type], + gen_key => receive_values.data, + ) + end + end + end + aggregate_solution = OrderedDict() + if mpi.root + for first_stage_var in FIRST_STAGE_VARS + aggregate_solution[first_stage_var] = OrderedDict() + for gen_key in keys(sp_solution[first_stage_var]) + aggregate_solution[first_stage_var][gen_key] = + sp_solution[first_stage_var][gen_key] + end + end + for i in 1:mpi.nprocs + push!(aggregate_solution, "s$i" => OrderedDict()) + for (solution_key, solution_dict) in gather_solution + push!(aggregate_solution["s$i"], solution_key => OrderedDict()) + if solution_key !== "Spinning reserve (MW)" + for (gen_bus_key, values) in solution_dict + aggregate_solution["s$i"][solution_key][gen_bus_key] = + gather_solution[solution_key][gen_bus_key][(i-1)*T+1:i*T] + end + else + for (reserve_name, reserve_dict) in solution_dict + push!( + aggregate_solution["s$i"][solution_key], + reserve_name => OrderedDict(), + ) + for (gen_key, values) in reserve_dict + aggregate_solution["s$i"][solution_key][reserve_name][gen_key] = + gather_solution[solution_key][reserve_name][gen_key][(i-1)*T+1:i*T] + end + end + end + end + end + end + return aggregate_solution +end diff --git a/src/solution/methods/ProgressiveHedging/structs.jl b/src/solution/methods/ProgressiveHedging/structs.jl new file mode 100644 index 0000000..b339813 --- /dev/null +++ b/src/solution/methods/ProgressiveHedging/structs.jl @@ -0,0 +1,130 @@ +# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +module ProgressiveHedging +using JuMP, MPI, TimerOutputs +import ..SolutionMethod + +mutable struct TerminationCriteria + max_iterations::Int + max_time::Float64 + min_feasibility::Float64 + min_improvement::Float64 + min_iterations::Int + + function TerminationCriteria(; + max_iterations::Int = 1000, + max_time::Float64 = 14400.0, + min_feasibility::Float64 = 1e-3, + min_improvement::Float64 = 1e-3, + min_iterations::Int = 2, + ) + return new( + max_iterations, + max_time, + min_feasibility, + min_improvement, + min_iterations, + ) + end +end + +Base.@kwdef mutable struct IterationInfo + it_num::Int + sp_consensus_vals::Array{Float64,1} + global_consensus_vals::Array{Float64,1} + sp_obj::Float64 + global_obj::Float64 + it_time::Float64 + total_elapsed_time::Float64 + global_residual::Float64 + global_infeas::Float64 +end + +mutable struct Method <: SolutionMethod + consensus_vars::Union{Array{VariableRef,1},Nothing} + weights::Union{Array{Float64,1},Nothing} + initial_global_consensus_vals::Union{Array{Float64,1},Nothing} + num_of_threads::Int + ρ::Float64 + λ_default::Float64 + print_interval::Int + termination_criteria::TerminationCriteria + + function Method(; + consensus_vars::Union{Array{VariableRef,1},Nothing} = nothing, + weights::Union{Array{Float64,1},Nothing} = nothing, + initial_global_consensus_vals::Union{Array{Float64,1},Nothing} = nothing, + num_of_threads::Int = 1, + ρ::Float64 = 1.0, + λ_default::Float64 = 0.0, + print_interval::Int = 1, + termination_criteria::TerminationCriteria = TerminationCriteria(), + ) + return new( + consensus_vars, + weights, + initial_global_consensus_vals, + num_of_threads, + ρ, + λ_default, + print_interval, + termination_criteria, + ) + end +end + +struct FinalResult + obj::Float64 + vals::Any + infeasibility::Float64 + total_iteration_num::Int + wallclock_time::Float64 +end + +struct SpResult + obj::Float64 + vals::Array{Float64,1} +end + +Base.@kwdef mutable struct SubProblem + mip::JuMP.Model + obj::AffExpr + consensus_vars::Array{VariableRef,1} + weights::Array{Float64,1} +end + +Base.@kwdef struct SpSolution + obj::Float64 + consensus_vals::Array{Float64,1} + residuals::Array{Float64,1} +end + +Base.@kwdef mutable struct SpParams + ρ::Float64 + λ::Array{Float64,1} + global_consensus_vals::Array{Float64,1} +end + +struct MpiInfo + comm::Any + rank::Int + root::Bool + nprocs::Int + + function MpiInfo(comm) + rank = MPI.Comm_rank(comm) + 1 + is_root = (rank == 1) + nprocs = MPI.Comm_size(comm) + return new(comm, rank, is_root, nprocs) + end +end + +Base.@kwdef struct Callbacks + before_solve_subproblem::Any + after_solve_subproblem::Any + after_iteration::Any +end + +end