Merge remote-tracking branch 'rael/wecc' into hydro_scenario

Author: staadecker

docs/Performance.md

@@ -0,0 +1,74 @@
+# Performance
+
+Memory use and solve time are two important factors that we try to keep to a minimum in our models. There are multiple
+things one can do to improve performance.
+
+## Solving methods
+
+By far the biggest factor that impacts performance is the method used by Gurobi. The fastest method is barrier solve
+without crossover (use `--recommended-fast`)
+however this method often returns a suboptimal solution. The next fastest is barrier solve followed by crossover and
+simplex (use `--recommended`) which almost always works. In some cases barrier solve encounters numerical issues (
+see [`Numerical Issues.md`](./Numerical%20Issues.md))
+in which case the slower Simplex method must be used (`--solver-options-string method=1`).
+
+## Solver interface
+
+Solver interfaces are how Pyomo communicates with Gurobi (or any solver).
+
+There are two solver interfaces that you should know about: `gurobi` and `gurobi_direct`.
+
+- When using `gurobi`, Pyomo will write the entire model to a temporary text file and then start a *separate Gurobi
+  process* that will read the file, solve the model and write the results to another temporary text file. Once Gurobi
+  finishes writing the results Pyomo will read the results text file and load the results back into the Python program
+  before running post_solve (e.g. generate csv files, create graphs, etc). Note that these temporary text files are
+  stored in `/tmp` but if you use `--recommended-debug` Pyomo and Gurobi will instead use a `temp` folder in your model.
+
+- `gurobi_direct` uses Gurobi's Python library to create and solve the model directly in Python without the use of
+  intermediate text files.
+
+In theory `gurobi_direct` should be faster and more efficient however in practice we find that that's not the case. As
+such we recommend using `gurobi` and all our defaults do so. If someone has the time they could profile `gurobi_direct`
+to improve performance at which point we could make `gurobi_direct` the default (and enable `--save-warm-start` by default, see below).
+
+The `gurobi` interface has the added advantage of separating Gurobi and Pyomo into separate threads. This means that
+while Gurobi is solving and Pyomo is idle, the operating system can automatically move Pyomo's memory usage
+to [virtual memory](https://serverfault.com/questions/48486/what-is-swap-memory)
+which will free up more memory for Gurobi.
+
+## Warm starting
+
+Warm starting is the act of using a solution from a previous similar model to start the solver closer to your expected
+solution. Theoretically this can help performance however in practice there are several limitations. For this section, *
+previous solution* refers to the results from an already solved model that you are using to warm start the solver. *
+Current solution* refers to the solution you are trying to find while using the warm start feature.
+
+- To warm start a model use `switch solve --warm-start <path_to_previous_solution>`.
+
+- Warm starting only works if the previous solution does not break any constraints of the current solution. This usually
+  only happens if a) the model has the exact same set of variables b)
+  the previous solution was "harder" (e.g. it had more constraints to satisfy).
+
+- Warm starting always uses the slower Simplex method. This means unless you expect the previous solution and current
+  solution to be very similar, it may be faster to solve without warm start using the barrier method.
+
+- If your previous solution didn't use crossover (e.g. you used `--recommended-fast`) then warm starting will be even
+  slower since the solver will need to first run crossover before warm starting.
+
+- Our implementation of warm starting only works if your previous solution has an `outputs/warm_start.pickle`
+  file. This file is only generated when you use `--save-warm-start`.
+
+- `--save-warm-start` and `--warm-start` both use an extension of the `gurobi_direct` solver interface which is
+  generally slower than the `gurobi` solver interface (see section above).
+  
+## Tools for improving performance
+
+- [Memory profiler](https://pypi.org/project/memory-profiler/) for generating plots of the memory
+use over time. Use `mprof run --interval 60 --multiprocess switch solve ...` and once solving is done
+  run `mprof plot -o profile.png` to make the plot.
+  
+- [Fil Profiler](https://pypi.org/project/filprofiler/) is an amazing tool for seeing which parts of the code are
+using up memory during peak memory usage.
+  
+- Using `switch_model.utilities.StepTimer` to measure how long certain code blocks take to run. See examples
+throughout the code.

switch_model/balancing/load_zones.py

@@ -286,9 +286,9 @@ def get_component_per_year(m, z, p, component):
 @graph(
     "energy_balance_duals",
     title="Energy balance duals per period",
-    note="Note: Outliers and zero-valued duals are ignored."
+    note="Note: Outliers and zero-valued duals are ignored.",
 )
-def graph(tools):
+def graph_energy_balance(tools):
     load_balance = tools.get_dataframe("load_balance.csv")
     load_balance = tools.transform.timestamp(load_balance)
     load_balance["energy_balance_duals"] = (
@@ -303,22 +303,50 @@ def graph(tools):
     # Don't include the zero-valued duals
     load_balance = load_balance.replace(0, tools.np.nan)
     if load_balance.count().sum() != 0:
-<<<<<<< HEAD
-        ax = tools.get_axes(
-            "energy_balance_duals",
-            title="Energy balance duals per period",
-            note="Note: Outliers and zero-valued duals are ignored.",
-        )
         load_balance.plot.box(
-            ax=ax,
+            ax=tools.get_axes(),
             xlabel="Period",
             ylabel="Energy balance duals (cents/kWh)",
             showfliers=False,
-=======
-        load_balance.plot.box(
-            ax=tools.get_axes(),
-            xlabel='Period',
-            ylabel='Energy balance duals (cents/kWh)',
-            showfliers=False
->>>>>>> b3590fdb (Redesign graphing API)
         )
+
+
+@graph("daily_demand", title="Total daily demand", supports_multi_scenario=True)
+def demand(tools):
+    df = tools.get_dataframe("loads.csv", from_inputs=True, drop_scenario_info=False)
+    df = df.groupby(["TIMEPOINT", "scenario_name"], as_index=False).sum()
+    df = tools.transform.timestamp(df, key_col="TIMEPOINT", use_timepoint=True)
+    df = df.groupby(
+        ["season", "hour", "scenario_name", "time_row"], as_index=False
+    ).mean()
+    df["zone_demand_mw"] /= 1e3
+    pn = tools.pn
+
+    plot = (
+        pn.ggplot(df)
+        + pn.geom_line(pn.aes(x="hour", y="zone_demand_mw", color="scenario_name"))
+        + pn.facet_grid("time_row ~ season")
+        + pn.labs(x="Hour (PST)", y="Demand (GW)", color="Scenario")
+    )
+    tools.save_figure(plot.draw())
+
+
+@graph("demand", title="Total demand", supports_multi_scenario=True)
+def yearly_demand(tools):
+    df = tools.get_dataframe("loads.csv", from_inputs=True, drop_scenario_info=False)
+    df = df.groupby(["TIMEPOINT", "scenario_name"], as_index=False).sum()
+    df = tools.transform.timestamp(df, key_col="TIMEPOINT", use_timepoint=True)
+    df["zone_demand_mw"] *= df["tp_duration"] / 1e3
+    df["day"] = df["datetime"].dt.day_of_year
+    df = df.groupby(["day", "scenario_name", "time_row"], as_index=False)[
+        "zone_demand_mw"
+    ].sum()
+    pn = tools.pn
+
+    plot = (
+        pn.ggplot(df)
+        + pn.geom_line(pn.aes(x="day", y="zone_demand_mw", color="scenario_name"))
+        + pn.facet_grid("time_row ~ .")
+        + pn.labs(x="Day of Year", y="Demand (GW)", color="Scenario")
+    )
+    tools.save_figure(plot.draw())

switch_model/generators/core/build.py

@@ -720,22 +720,7 @@ def post_solve(m, outdir):
     )
 
 
-<<<<<<< HEAD
-def graph(tools):
-    graph_capacity(tools)
-    graph_buildout_per_tech(tools)
-
-
-def compare(tools):
-    graph_buildout(tools)
-
-
-=======
-@graph(
-    "generation_capacity_per_period",
-    title="Online Generation Capacity Per Period"
-)
->>>>>>> b3590fdb (Redesign graphing API)
+@graph("generation_capacity_per_period", title="Online Generation Capacity Per Period")
 def graph_capacity(tools):
     # Load gen_cap.csv
     gen_cap = tools.get_dataframe("gen_cap.csv")
@@ -768,33 +753,25 @@ def graph_capacity(tools):
 
     # Plot
     # Get a new set of axis to create a breakdown of the generation capacity
-<<<<<<< HEAD
-    ax = tools.get_axes(
-        out="generation_capacity_per_period",
-        title="Online generating capacity by period",
-    )
     capacity_df.plot(
         kind="bar",
-        ax=ax,
-=======
-    capacity_df.plot(
-        kind='bar',
         ax=tools.get_axes(),
->>>>>>> b3590fdb (Redesign graphing API)
         stacked=True,
         ylabel="Capacity Online (GW)",
         xlabel="Period",
         color=tools.get_colors(len(capacity_df.index)),
     )
 
+    tools.bar_label()
+
 
 @graph(
     "buildout_gen_per_period",
     title="Built Capacity per Period",
-    supports_multi_scenario=True
+    supports_multi_scenario=True,
 )
 def graph_buildout(tools):
-    build_gen = tools.get_dataframe("BuildGen.csv")
+    build_gen = tools.get_dataframe("BuildGen.csv", dtype={"GEN_BLD_YRS_1": str})
     build_gen = build_gen.rename(
         {
             "GEN_BLD_YRS_1": "GENERATION_PROJECT",
@@ -804,13 +781,7 @@ def graph_buildout(tools):
         axis=1,
     )
     build_gen = tools.transform.build_year(build_gen)
-<<<<<<< HEAD
-    gen = tools.get_dataframe(
-        "generation_projects_info", from_inputs=True, all_scenarios=True
-    )
-=======
     gen = tools.get_dataframe("generation_projects_info", from_inputs=True)
->>>>>>> b3590fdb (Redesign graphing API)
     gen = tools.transform.gen_type(gen)
     gen = gen[["GENERATION_PROJECT", "gen_type", "scenario_name"]]
     build_gen = build_gen.merge(
@@ -845,33 +816,22 @@ def graph_buildout(tools):
 
     # Plot
     # Get a new set of axis to create a breakdown of the generation capacity
-<<<<<<< HEAD
-    ax = tools.get_axes(out="buildout_per_period", title="Built capacity per period")
     build_gen.plot(
         kind="bar",
-        ax=ax,
-=======
-    build_gen.plot(
-        kind='bar',
         ax=tools.get_axes(),
->>>>>>> b3590fdb (Redesign graphing API)
         stacked=True,
         ylabel="Capacity Online (GW)",
         xlabel="Period",
         color=tools.get_colors(len(build_gen.index)),
     )
-<<<<<<< HEAD
-
-=======
->>>>>>> b3590fdb (Redesign graphing API)
 
 
 @graph(
     "gen_buildout_per_tech_period",
     title="Buildout relative to max allowed for period",
     note="\nNote 1: This graph excludes predetermined buildout and projects that have no capacity limit."
-         "\nTechnologies that contain projects with no capacity limit are marked by a * and their graphs may"
-         "be misleading."
+    "\nTechnologies that contain projects with no capacity limit are marked by a * and their graphs may"
+    "be misleading.",
 )
 def graph_buildout_per_tech(tools):
     # Load gen_cap.csv
@@ -920,32 +880,16 @@ def graph_buildout_per_tech(tools):
     # Set the name of the legend.
     df = df.rename_axis("Type", axis="columns")
     # Add a * to tech
-<<<<<<< HEAD
     df = df.rename(
         lambda c: f"{c}*" if c in unlimited_gen_types.values else c, axis="columns"
     )
-    # Get axes to graph on
-    ax = tools.get_axes(
-        out="gen_buildout_per_tech_no_pred",
-        title="Buildout relative to max allowed for period",
-        note="\nNote 1: This graph excludes predetermined buildout and projects that have no capacity limit."
-        "\nTechnologies that contain projects with no capacity limit are marked by a * and their graphs may"
-        "be misleading.",
-    )
-=======
-    df = df.rename(lambda c: f"{c}*" if c in unlimited_gen_types.values else c, axis='columns')
->>>>>>> b3590fdb (Redesign graphing API)
     # Plot
     colors = tools.get_colors()
     if colors is not None:
         # Add the same colors but with a * to support our legend.
         colors.update({f"{k}*": v for k, v in colors.items()})
-<<<<<<< HEAD
-    df.plot(ax=ax, kind="line", color=colors, xlabel="Period", marker="x")
-=======
     ax = tools.get_axes()
-    df.plot(ax=ax, kind='line', color=colors, xlabel='Period', marker="x")
->>>>>>> b3590fdb (Redesign graphing API)
+    df.plot(ax=ax, kind="line", color=colors, xlabel="Period", marker="x")
     # Set the y-axis to use percent
     ax.yaxis.set_major_formatter(tools.mplt.ticker.PercentFormatter(1.0))
     # Horizontal line at 100%