Skip to content

Commit

Permalink
Updated solar hybridization example for EAGO v0.7.1
Browse files Browse the repository at this point in the history 
Note that in the lower_problem! and upper_problem! definitions, if the sense of the @NLobjective is changed to "Min", the "flcs" term in @NLobjective and the two opt._lower_objective_value/opt._upper_objective_value definitions also need to be negated. Base EAGO is assuming this is a minimization problem, so if the lower/upper problems return values consistent with a maximization problem, it will cause errors.
  • Loading branch information
@RXGottlieb
RXGottlieb committed Jul 5, 2022
1 parent 6024984 commit e78d2cb
Show file tree
Hide file tree
Showing 2 changed files with 111 additions and 84 deletions.
1 change: 1 addition & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
EAGO = "bb8be931-2a91-5aca-9f87-79e1cb69959a"
Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Expand Down
194 changes: 110 additions & 84 deletions notebooks/SolarHybridization/SolarHybridization.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
"metadata": {},
"outputs": [],
"source": [
"using JuMP, EAGO, Ipopt;"
"using JuMP, EAGO, Ipopt"
]
},
{
Expand All @@ -42,41 +42,44 @@
"metadata": {},
"outputs": [],
"source": [
"using CSV, DataFrames;"
"using CSV, DataFrames"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The code is organized as follows: \n",
"1. We import the solar resource data for the region we are concerned with. By default, this is the typical meteorological year (TMY) hourly data (8760 data points) which includes the geographic information (coordinates and timezone) as well as the direct normal irradiance (DNI). Together, we can model the process of concentrating the solar radiation and collecting it as heat in our system. \n",
"2. Once we have the resource data, we call \"PTCmodel.jl\" to simulate the specific thermal power of the solar concentrating technology (a parabolic trough in this case) in units of kW/m$^2$ given our geographic region and incident angle data. This function calls \"solarAngles.jl\" to calculate the angles of the direct solar radiation incident to the concentrator aperture with respect to each hour in the TMY data.\n",
"4. We define a function closure for the solar fraction calculated by \"iphProcessSmooth.jl\", which simulates the performance of the full solar concentrator and thermal storage system hybridized with the industrial process heat system. The solar fraction is calculated by Eq. 12 in the paper.\n",
"1. We import the solar resource data for the region we are concerned with. By default, this is the typical meteorological year (TMY) hourly data (8760 data points) which includes the geographic information (coordinates and timezone) as well as the direct normal irradiance (DNI). Combined, we can model the process of concentrating the solar radiation and collecting it as heat in our system. \n",
"2. Once we have the resource data, we call the `PTCmodel` function from \"PTCmodel.jl\" to simulate the specific thermal power of the solar concentrating technology (a parabolic trough in this case) in units of kW/m$^2$ given our geographic region and incident angle data. This function calls `solarAngles` from \"solarAngles.jl\" to calculate the angles of the direct solar radiation incident to the concentrator aperture with respect to each hour in the TMY data.\n",
"4. We define a function closure for the solar fraction calculated by `iphProcessSmooth` from \"iphProcessSmooth.jl\", which simulates the performance of the full solar concentrator and thermal storage system hybridized with the industrial process heat system. The solar fraction is calculated by Eq. 12 in the paper.\n",
"5. We define the custom upper- and lower-bounding optimization subproblems for the spatial branch-and-bound algorithm. The lower-bounding problem uses the convex capital model for its objective function (Eq. 16 in the paper), and the upper-bounding problem uses the nonconvex capital model as its objective function (Eq. 14 in the paper).\n",
"6. We set up the JuMP model with the EAGO optimizer (with custom bounding procedures) and we solve it."
]
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"include(\"smoothMinMaxAbs.jl\")\n",
"include(\"solarAngles.jl\")\n",
"include(\"PTCmodel.jl\")\n",
"include(\"iphProcessSmooth.jl\")\n",
"include(\"lifecycleCost.jl\")\n",
"# Step 1: read the data into a table and extract the appropriate data into a vector\n",
"include(\"smoothMinMaxAbs.jl\")\n",
"\n",
"# Step 1: Read the data into a table and extract the appropriate data into a vector\n",
"input_file = \"FirebaughTMY_Julia.csv\"\n",
"solData = CSV.File(input_file) |> DataFrame\n",
"yData = convert(Array{Float64,1},solData[:,7])\n",
"yData = convert(Array{Float64,1}, solData[:,7])\n",
"\n",
"# Step 2: get the specific thermal power potential for the region [kW/m^2]\n",
"q = PTCmodel(yData,-8,36.85,-120.46)\n",
"# Step 2: Get the specific thermal power potential for the region [kW/m^2]\n",
"q = PTCmodel(yData, -8, 36.85, -120.46)\n",
"\n",
"# Step 3: define the solar fraction closure\n",
"SolarFrac(xts,xa) = iphProcessSmooth(q,xts,xa);"
"# Step 3: Define the solar fraction closure\n",
"SolarFrac(xts,xa) = iphProcessSmooth(q, xts, xa);"
]
},
{
Expand All @@ -88,10 +91,11 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import EAGO: Optimizer, GlobalOptimizer\n",
"struct SolarExt <: EAGO.ExtensionType end"
]
},
Expand All @@ -104,48 +108,54 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Lower Problem Definition\n",
"import EAGO.lower_problem!\n",
"function lower_problem!(t::SolarExt,opt::EAGO.Optimizer)\n",
"function lower_problem!(t::SolarExt, opt::GlobalOptimizer)\n",
" \n",
" # Get active node\n",
" n = opt._current_node\n",
" \n",
" # Creates model, adds variables, and register nonlinear expressions\n",
" m = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer,tol=1.0e-3,print_level=0))\n",
" xL = n.lower_variable_bounds; xU = n.upper_variable_bounds\n",
" @variable(m, xL[i] <= x[i=1:2] <= xU[i])\n",
" # define convex lifecycle savings objective function cover\n",
" # Create model and add variables\n",
" mL = JuMP.Model(JuMP.optimizer_with_attributes(Ipopt.Optimizer,\n",
" \"tol\" => 1.0e-3,\n",
" \"print_level\" => 0))\n",
" xL = n.lower_variable_bounds\n",
" xU = n.upper_variable_bounds\n",
" @variable(mL, xL[i] <= x[i=1:2] <= xU[i])\n",
" \n",
" # Define convex lifecycle savings for use in the objective function\n",
" flcsC(xts,xa) = lifecycleCost(q, xL, xU, xts, xa, :convex)\n",
" JuMP.register(m, :flcsC, 2, flcsC, autodiff=true)\n",
" JuMP.register(m, :SolarFrac, 2, SolarFrac, autodiff=true)\n",
"\n",
" # Define nonlinear function\n",
" @NLobjective(m, Min, -flcsC(x[1], x[2]))\n",
" #@NLobjective(m, Max, flcsC(x[1], x[2]))\n",
" @NLconstraint(m, g1, SolarFrac(x[1], x[2]) >= 0.0)# declare constraints\n",
" JuMP.optimize!(m)\n",
" \n",
" # Register the convex lifecycle savings and solar fraction closure with our model\n",
" JuMP.register(mL, :flcsC, 2, flcsC, autodiff=true)\n",
" JuMP.register(mL, :SolarFrac, 2, SolarFrac, autodiff=true)\n",
"\n",
" # Get primal status, termination status, determine if a global solution was obtained\n",
" termination_status = JuMP.termination_status(m)\n",
" primal_status = JuMP.primal_status(m)\n",
" feasible_flag = EAGO.is_feasible_solution(termination_status, primal_status)\n",
" # Define the nonlinear objective and constraint\n",
" @NLobjective(mL, Max, flcsC(x[1], x[2]))\n",
" @NLconstraint(mL, g1, SolarFrac(x[1], x[2]) >= 0.0) # Declare constraints\n",
" \n",
" # Optimize the model\n",
" JuMP.optimize!(mL)\n",
"\n",
" # Interpret status codes for branch-and-bound\n",
" if feasible_flag\n",
" opt._lower_objective_value = 1.0*(JuMP.objective_value(m)-1e-4)\n",
" tstatus = MOI.get(mL, MOI.TerminationStatus())\n",
" pstatus = MOI.get(mL, MOI.PrimalStatus())\n",
" \n",
" if EAGO.local_problem_status(tstatus, pstatus) == EAGO.LRS_FEASIBLE\n",
" opt._lower_objective_value = -1.0*(JuMP.objective_value(mL)-1e-4) # Multiplied by -1 because EAGO expects \"Min\"\n",
" opt._lower_solution = JuMP.value.(x)\n",
" opt._lower_feasibility = true\n",
" opt._cut_add_flag = false\n",
" else\n",
" opt._lower_feasibility = false\n",
" opt._lower_objective_value = -Inf\n",
" opt._lower_objective_value = Inf\n",
" opt._cut_add_flag = false\n",
" end\n",
" \n",
" return\n",
"end;"
]
Expand All @@ -159,45 +169,48 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Upper Problem Definition\n",
"import EAGO.upper_problem!\n",
"function upper_problem!(t::SolarExt,opt::EAGO.Optimizer)\n",
"function upper_problem!(t::SolarExt,opt::GlobalOptimizer)\n",
" \n",
" # Get active node\n",
" n = opt._current_node\n",
" \n",
" # Creates model, adds variables, and register nonlinear expressions\n",
" m = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer,tol=1.0e-3,print_level=0))\n",
" xL = n.lower_variable_bounds; xU = n.upper_variable_bounds\n",
" @variable(m, xL[i] <= x[i=1:2] <= xU[i])\n",
" # define lifecycle savings objective function cover\n",
" # Create model and add variables\n",
" mU = JuMP.Model(JuMP.optimizer_with_attributes(Ipopt.Optimizer,\n",
" \"tol\" => 1.0e-3,\n",
" \"print_level\" => 0))\n",
" xL = n.lower_variable_bounds\n",
" xU = n.upper_variable_bounds\n",
" @variable(mU, xL[i] <= x[i=1:2] <= xU[i])\n",
" \n",
" # Define nonconvex lifecycle savings for use in the objective function\n",
" flcs(xts,xa) = lifecycleCost(q, xL, xU, xts, xa, :nonconvex)\n",
" JuMP.register(m, :flcs, 2, flcs, autodiff=true)\n",
" JuMP.register(m, :SolarFrac, 2, SolarFrac, autodiff=true)\n",
" \n",
" # Register the convex lifecycle savings and solar fraction closure with our model\n",
" JuMP.register(mU, :flcs, 2, flcs, autodiff=true)\n",
" JuMP.register(mU, :SolarFrac, 2, SolarFrac, autodiff=true)\n",
"\n",
" # Define nonlinear function\n",
" @NLobjective(m, Min, -flcs(x[1], x[2]))\n",
" #@NLobjective(m, Max, flcs(x[1], x[2]))\n",
" @NLconstraint(m, g1, SolarFrac(x[1], x[2]) >= 0.0)# declare constraints\n",
" JuMP.optimize!(m)\n",
" @NLobjective(mU, Max, flcs(x[1], x[2]))\n",
" @NLconstraint(mU, g1, SolarFrac(x[1], x[2]) >= 0.0)# declare constraints\n",
" JuMP.optimize!(mU)\n",
"\n",
" # Get primal status, termination status, determine if a global solution was obtained\n",
" termination_status = JuMP.termination_status(m)\n",
" primal_status = JuMP.primal_status(m)\n",
" feasible_flag = EAGO.is_feasible_solution(termination_status, primal_status)\n",
" # Interpret status codes for branch-and-bound\n",
" tstatus = MOI.get(mU, MOI.TerminationStatus())\n",
" pstatus = MOI.get(mU, MOI.PrimalStatus())\n",
"\n",
" # Interpret status codes for branch and bound\n",
" if feasible_flag\n",
" opt._upper_objective_value = 1.0*JuMP.objective_value(m)\n",
" if EAGO.local_problem_status(tstatus, pstatus) == EAGO.LRS_FEASIBLE\n",
" opt._upper_objective_value = -1.0*JuMP.objective_value(mU) # Multiplied by -1 because EAGO expects \"Min\"\n",
" opt._upper_solution = JuMP.value.(x)\n",
" opt._upper_feasibility = true\n",
" else\n",
" opt._upper_feasibility = false\n",
" opt._upper_objective_value = -Inf\n",
" opt._upper_objective_value = Inf\n",
" opt._cut_add_flag = false\n",
" end\n",
" return\n",
Expand All @@ -213,42 +226,48 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"import EAGO: preprocess!, postprocess!, cut_condition\n",
"function EAGO.preprocess!(t::SolarExt, x::Optimizer)\n",
"function EAGO.preprocess!(t::SolarExt, x::GlobalOptimizer)\n",
" x._preprocess_feasibility = true\n",
" return\n",
"end\n",
"function EAGO.postprocess!(t::SolarExt, x::Optimizer)\n",
"function EAGO.postprocess!(t::SolarExt, x::GlobalOptimizer)\n",
" x._postprocess_feasibility = true\n",
" return\n",
"end\n",
"EAGO.cut_condition(t::SolarExt, x::Optimizer) = false"
"EAGO.cut_condition(t::SolarExt, x::GlobalOptimizer) = false"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we define the JuMP model and the variables (and bounds). We must specify that our new custom extension $\\texttt{SolarExt}$ to EAGO's default routines should be used and that we will be branching on both of our decision variables. The latter is required for custom routines since no expressions will be provided to the EAGO optimizer and therefore it cannot infer which variables should be branched on. "
"Now, we define the JuMP model and the variables (and bounds). We must specify that our new custom extension $\\texttt{SolarExt}$ of EAGO's default routines should be used and that we will be branching on both of our decision variables. The latter is required for custom routines since no expressions will be provided to the EAGO optimizer and therefore it cannot infer which variables should be branched on. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"m = JuMP.Model(with_optimizer(EAGO.Optimizer,relative_tolerance=1e-2,verbosity=1,\n",
" output_iterations=1, branch_variable=Bool[true; true],\n",
" ext_type=SolarExt()))\n",
"factory = () -> EAGO.Optimizer(SubSolvers(; t=SolarExt()))\n",
"m = JuMP.Model(optimizer_with_attributes(factory,\n",
" \"relative_tolerance\" => 1e-2,\n",
" \"verbosity\" => 1,\n",
" \"output_iterations\" => 1, \n",
" \"branch_variable\" => Bool[true; true],\n",
" ))\n",
"\n",
"x_L = [0.001,000.01]\n",
"x_U = [16.0,60000.0]\n",
"@variable(m, x_L[i] <= x[i=1:2] <= x_U[i]);"
"x_L = [0.001, 0.01]\n",
"x_U = [16.0, 60000.0]\n",
"@variable(m, x_L[i] <= x[i=1:2] <= x_U[i])\n",
"JuMP.register(m, :SolarFrac, 2, SolarFrac, autodiff=true)\n",
"@NLconstraint(m, g1, SolarFrac(x[1], x[2]) >= 0.0);"
]
},
{
Expand All @@ -260,37 +279,44 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"******************************************************************************\n",
"This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
" Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
" For more information visit https://github.com/coin-or/Ipopt\n",
"******************************************************************************\n",
"\n",
"-----------------------------------------------------------------------------------------------------------------------------\n",
"| Iteration # | Nodes | Lower Bound | Upper Bound | Gap | Ratio | Time | Time Left |\n",
"-----------------------------------------------------------------------------------------------------------------------------\n",
"| 1 | 2 | --0.753E7 | 0.308E3 | 0.753E7 | 0.100E1 | 0.647E0 | 0.999E3 |\n",
"| 2 | 3 | --0.753E7 | 0.308E3 | 0.753E7 | 0.100E1 | 0.142E2 | 0.986E3 |\n",
"| 3 | 4 | --0.749E7 | 0.308E3 | 0.749E7 | 0.100E1 | 0.279E2 | 0.972E3 |\n",
"| 4 | 5 | --0.749E7 | 0.308E3 | 0.749E7 | 0.100E1 | 0.408E2 | 0.959E3 |\n",
"| 5 | 2 | --0.734E7 | --0.732E7 | 0.265E5 | 0.361E-2 | 0.420E2 | 0.958E3 |\n",
"| 1 | 2 | -7.531E+06 | 3.084E+02 | 7.531E+06 | 1.000E+00 | 9.286E+00 | 3.591E+03 |\n",
"| 2 | 3 | -7.531E+06 | 3.084E+02 | 7.531E+06 | 1.000E+00 | 2.178E+01 | 3.578E+03 |\n",
"| 3 | 4 | -7.492E+06 | 3.084E+02 | 7.492E+06 | 1.000E+00 | 3.428E+01 | 3.566E+03 |\n",
"| 4 | 5 | -7.492E+06 | 3.084E+02 | 7.492E+06 | 1.000E+00 | 4.621E+01 | 3.554E+03 |\n",
"| 5 | 6 | -7.344E+06 | -7.318E+06 | 2.654E+04 | 3.614E-03 | 4.728E+01 | 3.553E+03 |\n",
" \n",
"Relative Tolerance Achieved\n",
"First Solution Found at Node 9\n",
"LBD = -7.344210617924415e6\n",
"UBD = -7.3176663129719095e6\n",
"LBD = -7.344210617732814e6\n",
"UBD = -7.3176663129717745e6\n",
"Solution is :\n",
" X[1] = 11.722361246205823\n",
" X[2] = 43615.19059009129\n",
" X[1] = 11.72236127830526\n",
" X[2] = 43615.19078439399\n",
" \n",
" 41.952241 seconds (5.48 M allocations: 60.761 GiB, 4.64% gc time)\n",
"xts* = 11.722361246205823 xa* = 43615.19059009129 f* = 7.3176663129719095e6 SF* = 0.6980196900769012\n",
"xts* = 11.72236127830526 xa* = 43615.19078439399 f* = 7.3176663129717745e6 SF* = 0.698019691472141\n",
"Algorithm terminated with a status of OPTIMAL and a result code of FEASIBLE_POINT\n"
]
}
],
"source": [
"@time JuMP.optimize!(m)\n",
"optimize!(m)\n",
"\n",
"println(\"xts* = \", JuMP.value(x[1]), \" xa* = \",\n",
" JuMP.value(x[2]),\" f* = \",-1.0*JuMP.objective_value(m),\" SF* = \",\n",
Expand All @@ -314,15 +340,15 @@
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.6.1",
"display_name": "Julia 1.7.3",
"language": "julia",
"name": "julia-1.6"
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.1"
"version": "1.7.3"
}
},
"nbformat": 4,
Expand Down

0 comments on commit e78d2cb

Please sign in to comment.