diff --git a/doc/theory_manual/ComponentCharacterization/HERON_Cashflows.md b/doc/theory_manual/ComponentCharacterization/HERON_Cashflows.md
new file mode 100644
index 00000000..dc27c4a1
--- /dev/null
+++ b/doc/theory_manual/ComponentCharacterization/HERON_Cashflows.md
@@ -0,0 +1,192 @@
+
+
+![HERON_cfs](../diagrams/HERON_cfs.png)
+
+## ***Cashflow Definition***
+
+Cash flows are defined to map resource-dependent component activities, capacities and market uncertainty to economic values. The engine for calculating economic metrics is the Tool for Economic Analysis (TEAL) which is a plugin in the FORCE toolset (more specifically, a plugin for RAVEN). Cash flows are defined by a universal formula:
+
+$$ \hat{F} = \alpha \Big(\frac{d}{d^\prime}\Big)^\chi$$
+
+where the individual terms are:
+
+- $\alpha$: representative price or cost
+- $d$: driver of the cash flow
+- $d^\prime$: reference driver associated with the cost, if needed (usually 1)
+- $\chi$: scaling factor representing economies ($0\leq\chi\leq1$) or diseconomies of scale ($\chi>1$)
+
+A given cash flow is defined by a set of $(\alpha, d, d^\prime, \chi)$ values. Cash flows with $\chi\neq1$ typically introduce nonlinearities to any optimization problem but are still solvable by some open source solvers like IPOPT. A lot of cash flows are represented simply as $(\alpha, d, 1, 1)$.
+
+The individual cash flow parameters $(\alpha, d, d^\prime, \chi)$ can be defined by static values or their value can be defined by other variables within the bi-level optimization. For example, a cash flow driver could be the capacity of a certain HERON Component (e.g., a plant or generator) or the dispatch activity determined at some time step (e.g., hourly).
+
+## **Additional Cash Flow Modifications**
+
+Cash flows can be modified based on tax rates and inflation rates per year of the project simulation. Tax and inflation can be applied on a per cash flow basis. Additionally a discount rate or weighted average cost of capital (WACC) can be introduced to determine the present value of future cash flows. Final cash flows can be summarized with the following formula for a given year $y$:
+
+$$ F_y = \frac{1-\tau}{(1+\nu)^y (1+r)^y}\hat{F} $$
+
+or, after some simplification,
+
+$$ F_y = \frac{1-\tau}{(1+\nu+r+\nu r)^y}\hat{F} $$
+
+where
+
+- $\tau$: tax rate
+- $\nu$: inflation rate
+- $r$: discount rate or WACC
+
+SInce tax and inflation is toggle-able for each defined cash flow in HERON, we can define indicator functions:
+
+$$\tau_{p}^F = 1-\delta_{p}^{\tau}\tau$$
+
+where
+
+$$ \delta_{p}^{\tau} =
+\begin{cases}
+1\hspace{0.5cm} \text{if tax is applied for component } p \text{ and cash flow } F \text{ and no depreciation}\\
+0\hspace{0.5cm} \text{elsewhere }
+\end{cases}$$
+
+and the inflation is applied similarly as
+
+$$\nu_{p,y}^F = (1+ \delta_{p}^{\nu}\nu)^{-y} = \frac{1}{(1+ \delta_{p}^{\nu}\nu)^y}.$$
+
+The user can specify whether each individual cashflow for each component are either taxable and/or should be inflated. With tax and inflation, we can introduce another shorthand with
+
+$$\lambda_{p,y}^F = \tau_{p}^F \nu_{p,y}^F$$
+
+which defaults to 1 in the case of no taxes or inflation. The final definition of cash flow is written as:
+$$ F_y = \frac{\lambda_{p,y}^F}{(1+r)^y}\hat{F} $$
+
+These get applied to the summed yearly cashflows. Note that cashflows are typically a function of dispatch or capacity in HERON:
+
+$$F_y = f(\mathbf{c}, {}^s\mathbf{D}^\star)$$
+
+## **Types of Cash Flows**
+
+Three distinct groups of cash flows are defined in TEAL:
+
+1. One-Time (and Depreciation)
+2. Recurring Yearly
+3. Recurring Hourly
+
+
+### *One-Time*
+One-time cash flows (e.g., capital expenditures or CAPEX) are applied once per lifetime of the component. Let's call the specific One-Time cashflow for Capital Expenditures, or CAPEX, $A$:
+
+$$ \hat{F}^{one-time}_p \equiv A_p =  \alpha_p^A \left(\frac{c_p}{d^{A\prime}_p} \right)^{\chi^A_p}$$
+
+Note that if component lifetimes are shorter than the total project simulation, these capital expenditures are repeated to simulate reconstruction costs. Reconstruction is assumed to happen instantaneously. We can use an indicator function to demonstrate this. Each component $p$ has a lifetime $L_p$. We define a set of integers dependent on whether the component lifetime is less than the project lifetime $Y$:
+
+$$\mathbb{N}_p \in \mathbb{Z}^+$$
+
+and
+
+$$\mathbb{N}_p =
+\begin{align*}
+\begin{cases}
+\{0\}\hspace{0.5cm} &\text{if } L_p > Y\\
+\{0, 1\}\hspace{0.5cm} &\text{if } L_p = Y\\
+\{0, 1, \ldots, \big\lfloor \frac{Y}{L_p} \big\rfloor \}\hspace{0.5cm} &\text{if } L_p < Y\\
+\end{cases} \end{align*}$$
+
+We define a function $\epsilon$ of the given year $y$ and component $p$ as
+
+$$ \epsilon_{p,y}^A =
+\begin{cases}
+1\hspace{0.5cm} \text{if } y = nL_p \ \forall \ n \in \mathbb{N}_p\\
+0\hspace{0.5cm} \text{elsewhere }
+\end{cases}$$
+
+which is equal to 1 when the year $y$ in consideration is a multiple of the component lifetime (starting with and including year 0). Then we can write the full CAPEX cost as a function not just of component but of year $y$:
+
+$$ A_{p,y} =  \epsilon_{p,y}^A \alpha_p^A \left(\frac{c_p}{c^{A\prime}_p} \right)^{\chi^A_p}.$$
+
+We can also differentiate between a positive and negative one-time cash flow:
+
+$$ A_{p,y}^+,  A_{p,y}^-$$
+
+#### *Depreciation*
+Users have the option to apply depreciation to the Cash Flows (only for capital expenditures, really). The value of an asset deteriorates over time, so (per the IRS of the United States) there is a system to recover the depreciation of the asset's value over a set recovery period through tax credits. TEAL offers yearly depreciation rates according to the MACRS (modified accelerated cost recovery system) for 3-, 5-, 7-, 10-, and 15-year periods. In reality it is recovery time + 1 yr due to convention: the asset on the books is shown to have been purchased in the middle of the year no matter when the purchase happened, so there is an extra payment a year after the final year in the recovery period.
+
+Say we select a recovery period $R_p$ for a component such that
+
+$$R_p \leq L_p \leq Y $$
+
+and
+
+$$R_p \in \{3,5,7,10,15\}. $$
+
+According to MACRS, we would have a depreciation rate $\beta_{p,y}$ per year. For example, for $R_p=3$ the rates would be
+
+$$\beta_{p,y}\Big|_{R_p=3} =
+\begin{cases}
+33.33\hspace{0.5cm} \ \ \text{if } y = nL_p + 1 \ \forall \ n \in \mathbb{N}_p\\
+44.45\hspace{0.5cm} \ \ \text{if } y = nL_p + 2 \ \forall \ n \in \mathbb{N}_p\\
+14.81\hspace{0.5cm} \ \ \text{if } y = nL_p + 3 \ \forall \ n \in \mathbb{N}_p\\
+7.41\hspace{0.7cm} \ \ \text{if } y = nL_p + 4\ \forall \ n \in \mathbb{N}_p\\
+0\hspace{0.75cm} \ \ \text{elsewhere }
+\end{cases}$$
+
+Note that no depreciation is applied during the construction year $y=0$. Users can alternatively apply custom depreciation rate schedules. To model a depreciation asset, TEAL applies two additional cash flows to the one-time CAPEX. First is a positive cash flow representing the tax credit that helps recuperate the depreciation. It is untaxed, but inflation *is* applied:
+
+$$B_{p,y}^+ = \beta_{p,y}\nu^B_{p,y}A_{p,0}.$$
+
+Note that within $A_{p,0}$ is the term $\epsilon_{p,y}^A$ which applies the capital cost ONLY on the first year of construction (and subsequent reconstructions). In TEAL, the driver of the cash flow is the entire first year cash flow. Since $R_p \leq L_p$ only a singular capital cost is applied per component lifetime. For a more complete definition, we would repeat the depreciation process again per lifetime (this is accounted for in the definition of $\beta$).
+
+A second, negative cash flow is also applied to represent the depreciation of the asset's value. This is similar to the first positive cash flow but it is taxed:
+
+$$B_{p,y}^- = \beta_{p,y}\lambda^B_{p,y}A_{p,0}.$$
+
+If we collect all like terms, the actual CAPEX cash flow will look like:
+
+$$\begin{align*} \sum_{j=\{A,B\},\ k=\{+,-\}} F_j^k = -\lambda_{p,y}A_{p,y}^A + B_{p,y}^+ - B_{p,y}^- &= -\lambda_{p,y}^A A_y + \beta_{p,y}\nu^B_{p,y}A_{p,0} - \beta_{p,y}\lambda^B_{p,y}A_{p,0} \\
+&= -\lambda^A_{p,y}A_y + \beta_{p,y}(\nu^B_{p,y} - \lambda^B_{p,y})A_{p,0} \\
+&= -\lambda^A_{p,y}A_y + \beta_{p,y}\nu^B_{p,y}(1 - \tau^B_{p,y})A_{p,0}
+\end{align*}$$
+
+If levelized cost of capital is selected, we should also collect the amortization and depreciation terms in the "multiplied" column and divide the remainder cashflow sum by those terms.
+
+
+### *Recurring Yearly*
+
+There are fixed yearly expenditures that are typically just indexed by component for all years after first "construction" year:
+
+```math
+\hat{F}_{p,y}^\text{yearly} \equiv \Gamma_{p,y>0} = \alpha_p^\Gamma \left(\frac{c_p}{c_p^{\Gamma\prime}} \right)^{\chi_p^\Gamma}
+```
+
+or
+
+$$\Gamma_{p,y} = \epsilon_y \alpha_p^{\Gamma} \left(\frac{c_p}{c_p^{\Gamma\prime}} \right)^{\chi_p^\Gamma} = \epsilon_y\alpha^{\Gamma}_p \left(\frac{c_p}{c_p^{\Gamma\prime}} \right)^{\chi_p^\Gamma}  $$
+
+where we define a parameter to zero-out terms in year 0
+
+$$\epsilon_{y} =
+\begin{cases}
+0\hspace{0.5cm} \ \text{if } y = 0 \ \forall \ y \in \mathbb{Y}\\
+1\hspace{0.5cm} \ \text{elsewhere }
+\end{cases}$$
+
+These could potentially have different costs per year (e.g., $\alpha^{\Gamma}_{p,y}$). But for now only indexed by component and driven by the component capacity.
+
+We can also differentiate between a positive and negative recurring yearly cash flow:
+
+$$ \Gamma_{p,y}^+,  \Gamma_{p,y}^-$$
+
+### *Recurring Hourly*
+
+The final type of cashflows, which is the most complex, is the hourly recurring cashflows. In reality, the "yearly" and "hourly" cash flows should be classified as pertaining to some macro step size (years) and a pivot parameter (at a smaller resolution, could be 15-min). When HERON runs dispatch optimization and collects all chosen dispatch strategies, it sends to TEAL a yearly sum of these cashflows. For our definition, we retain the subscript $y$ for the generic cashflow $F$ but write it as
+```math
+\hat{F}_{p,y}^{hourly} \equiv H_{p,y>0} = \sum\limits_{x=0}^X \sum\limits_{u=0}^U m_{y,u} \sum\limits_{t=0}^T \alpha_{p,x,y,u,t}^H \left( \frac{D_{p,x,y,u,t}}{D_{p,x,y,u,t}^{H \prime}} \right)^{\chi_{p,x,y,u,t}^H}
+```
+
+where $m_{y,u}$ is the multiplicity per cluster. This is an all-encompassing generic definition for an hourly cashflow. Note that the cashflow driver is no longer capacity, rather the dispatch (or production) of a given resource per component and time. Time is indexed further by year, cluster, and hour. Because it looks complex to explicitly specify all five indeces, we could instead use a notation
+
+$$i \equiv \{p,x,y,u,t\}$$
+
+and the $\epsilon_y$ parameter from the recurring yearly cash flow so that
+
+$$H_{p,y} = \epsilon_y \sum\limits_{x=0}^X \sum\limits_{u=0}^U m_{y,u} \sum\limits_{t=0}^T \alpha_{i}^H \left( \frac{D_{i}}{D_{i}^{H\prime}} \right)^{\chi_{i}^H}$$
+
+as a shorthand to demonstrate that we're indexing by all five indeces $\{p,x,y,u,t\}$. Note also that for every cluster, we are multiplying by its associated multiplicity: the amount of segments each cluster is meant to represent. This multiplicity may vary by year and cluster, but note that the number of clusters $U$ will not.
diff --git a/doc/theory_manual/ComponentCharacterization/HERON_Components.md b/doc/theory_manual/ComponentCharacterization/HERON_Components.md
new file mode 100644
index 00000000..fd3687ec
--- /dev/null
+++ b/doc/theory_manual/ComponentCharacterization/HERON_Components.md
@@ -0,0 +1,75 @@
+![HERON_components_resources](../diagrams/HERON_comps.png)
+
+## ***Component Capacities***
+We first define variables in our problem. Each component has the capability to perform three actions on a resource:
+  1) store,
+  2) demand,
+  3) produce or consume+produce.
+
+Each component can only perform one of those three actions. If it produces/consumes, then it can perform the action on multiple resources with a given transfer function that describes the conversion of a set $A$ of resources into another set $B$ of resources. For these given actions, the component will have a maximum value at which they can perform this action: these are capacities which, when collected for all components, is given as a vector
+$$\mathbf{c} = \big[ \ c_p : p \in \mathbb{P}^\prime\ \big]. $$
+
+Note that we are considering only the set of components that store or produce here, $\mathbb{P}^\prime \in \mathbb{P}$. This is closer to the definition of an IES; the difference is that within HERON we must include a resource sink for the resources to go to fulfill a demand or get sold. The demanding resources within HERON are typically markets or grids which are not typically a part of the IES but must be defined to conduct the simulation. We may be able to include the capacity of the market as an additional variable...
+
+## ***Component Minimums and Capacity Factors***
+There are optional inputs for all components to define the minimum level of production or resource activity. If not specified, the default is 0.
+$$\mathbf{m} = \big[ \ m_p : p \in \mathbb{P}^\prime\ \big]. $$
+
+Additionally, capacity factors can be used to modify the upper bound of production on an hourly time scale.
+$$\mathbf{f} = \big[ \ \big[\ \big[ \ \big[ f_{p,y,u,t} : t \in \mathbb{T} \ \big] : u \in \mathbb{U}_Y \ \big] : y \in \mathbb{Y}  \ \big]   : p \in \mathbb{P}^\prime\ \big]. $$
+
+## ***Transfer Functions***
+Only Producer components have the option to declare a transfer function which defines how a subset of resources is converted into another set. Three type of transfer functions are allowed:
+
+  1) Ratio/Linear
+  2) Polynomial
+  3) Function
+
+### *Ratio*
+The ratio transfer function mimics a chemical equation of the type
+$$3a + 7b \rightarrow 2c$$
+In this case, the production of resource $c$ is 3:2 with respect to resource $a$ and 7:2 with respect to resource $b$.
+
+### *Polynomial*
+Users can also define a polynomial relationship between the different resources. for example
+$$2a^2 + 3ab -b^2 + ba^3 = 0 $$
+
+### *Function*
+Custom Python methods are also allowed to be used to define the transfer relationship from one subset of resources to another.
+
+## ***Dispatch Activity***
+The actual action the component takes on resources per timestep is referred to as the dispatch activity and can be indexed similarly for each component by the resource being acted on as well as the time. We subdivide the time into years, then into clustered segments and then into a smaller timestep (typically hours per year).  The dispatch matrix is
+$${}^s\mathbf{D}^\star = \big[ \ \big[ \ \big[\ \big[ \ \big[ {}^sD^\star_{p,x,y,u,t} : t \in \mathbb{T} \ \big] : u \in \mathbb{U}_Y \ \big] : y \in \mathbb{Y}  \ \big] : x \in \mathbb{X} \ \big]: p \in \mathbb{P}^\prime \ \big]$$
+
+For the different component types, the dispatch activity is tracked via the following dispatch optimization variables:
+  1) Producer $\rightarrow$ production
+  2) Demand $\rightarrow$ production
+  3) Storage $\rightarrow$ [level, charge, discharge]
+
+The capacities are applied as a global upper bound on the dispatch when defined as flexible:
+
+$$ m_p \leq D_{p,x,y,u,t} \leq f_{p,y,u,t} c_p \ \forall \ \begin{cases} p \in \mathbb{P}^\prime \\ x\in \mathbb{X} \\ y\in \mathbb{Y} \\ u\in \mathbb{U}_Y \\ t\in \mathbb{T} \end{cases}  $$
+
+The minimum $m_p$ defaults to 0 unless specified; the capacity factor $f_{p,y,u,t}$ defaults to 1. In the default case, the bounds on dispatch are:
+
+$$ 0 \leq D_{p,x,y,u,t} \leq c_p \ \forall \ \begin{cases} p \in \mathbb{P}^\prime \\ x\in \mathbb{X} \\ y\in \mathbb{Y} \\ u\in \mathbb{U}{}_Y \\ t\in \mathbb{T} \end{cases}  $$
+
+Otherwise the components will just dispatch at a fixed level based on their capacity. We can also define the specific collection of dispatch actions taken within a specific scenario (i.e., for a given stochastic profile) as
+
+$${}^s\mathbf{D} \equiv \mathbf{D}(s)$$
+
+or perhaps more explicitly,
+
+```math
+{}^s\mathbf{D} \equiv \mathbf{D}({}^s\mathbf{W})
+```
+
+where we're not exactly indexing per scenario, it's more of a response of the dispatch to a given scenario $s$ and denote it as shown above.
+
+The optimal dispatch is defined as one which maximizes an inner objective function as shown below:
+
+```math
+{}^s\mathbf{D}^\star_{p,x,\hat{y},\hat{u},t} = \underset{D}{\text{arg max}} \ \sum_{P,X,T_{U}} {}^s\hat{F}_{hourly} \Big|_{\hat{u},\hat{y}}
+```
+
+Here, the objective is to maximize the sum of all *hourly* cashflows for all components, all resources, and all times per segment **for each** cluster/segment, **for each** year, **for each** scenario.
diff --git a/doc/theory_manual/HERON_Theory_Manual.md b/doc/theory_manual/HERON_Theory_Manual.md
new file mode 100644
index 00000000..8d861914
--- /dev/null
+++ b/doc/theory_manual/HERON_Theory_Manual.md
@@ -0,0 +1,29 @@
+# HERON Theory Manual
+
+Setting up a HERON run requires the following activities (taken from the HERON User Guide):
+
+1. Market Characterization
+   - Electricity and/or Commodity Markets
+   - Regulated vs. Deregulated Markets
+   - Price Taker vs. Price Maker modeling
+   - Grid-support functions (Ancilliary Services)
+   - Additional Tax Incentives
+2. Technology Identification
+3. Technology Characterization
+   - [(Technical) HERON Components](./ComponentCharacterization/HERON_Components.md)
+   - [(Economic) HERON CashFlows](./ComponentCharacterization/HERON_Cashflows.md)
+4. Time History Training
+   - [Time Index](./TimeHistoryTraining/HERON_TimeIndexing.md)
+5. Pre-analysis Calculations
+
+Within these activities, some additional guides are presented for HERON definitions.
+
+## HERON Workflows
+Different workflows:
+- `standard`: [(RAVEN-runs-RAVEN)](./Workflows/HERON_Standard_Workflow.md)
+  - most common, standard for a reason
+  - stochastic bi-level optimization scheme
+- `MOPED`:
+  - all-at-once solve (testing/development temporarily suspended)
+- `DISPATCHES`:
+  - all-at-once solve + integration with DISPATCHES algebraic models (development paused)
diff --git a/doc/theory_manual/TimeHistoryTraining/HERON_TimeIndexing.md b/doc/theory_manual/TimeHistoryTraining/HERON_TimeIndexing.md
new file mode 100644
index 00000000..7ed2a52d
--- /dev/null
+++ b/doc/theory_manual/TimeHistoryTraining/HERON_TimeIndexing.md
@@ -0,0 +1,43 @@
+
+![HERON_components_resources](../diagrams/HERON_time.png)
+
+
+### *Time Indexing*
+Consider a generator of stochastic profiles. We want $S$ profiles each indexed by $s \in \mathbb{S}$. Each profile consists of $Y$ years, each year with $T_Y$ total time steps. We can subdivide the year into $N$ segments each with a length of $T_U$. To do this, we first ensure that
+
+```math
+T_Y = N T_U
+```
+
+where $N$ is an integer number. One example of this is to have individual time steps of 1 hour, so that $T_Y = 8760$, and then divide the year into segments corresponding to days so that $T_U = 24$ and $N=365$ (ignoring leap years). One could also group instead by every 48 hours, every week, month, etc. as long as the number multiplies into $T_Y$ evenly.
+
+After training the synthetic history generator on each segment of each year, a clustering algorithm can determine how best to group the segments into $U$ unique groups each indexed by $u \in \mathbb{U}_Y$. If $\mathbb{N}_Y$ is the set of all segments per year, we would then have $U$ unique subsets of segments
+
+```math
+\mathbb{N}_u \subset \mathbb{N}_Y \ \ \forall u \in \mathbb{U}_Y
+```
+
+and
+
+```math\{\ \cup \ \mathbb{N}_u \ | u \in \mathbb{U}_Y \}=\mathbb{N}_Y
+```
+
+We can collect the number of segments each cluster represents as a multiplicity vector:
+
+```math$\mathbf{m} = [\ \ \text{len}(\mathbb{N}_u) \ | \ u \in \mathbb{U}_Y \ ]
+```
+
+The stochastic profiles which will be returned from the trained model will then have a shape of $Y \times U \times T$ and is given as
+
+```math
+{}^s\mathbf{W}
+```
+
+for a given evaluation $s$ of some stochastic profile generator. The value of the synthetic signal at some timestep $t$ of cluster $u$ during year $y$ is indexed as
+
+```math
+ {}^s\mathbf{W}_{y,u,t}
+ ```
+
+This stochastic profile could have an additional dimension if multiple signals are produced in one profile.
+
diff --git a/doc/theory_manual/Workflows/HERON_Standard_Workflow.md b/doc/theory_manual/Workflows/HERON_Standard_Workflow.md
new file mode 100644
index 00000000..37e341d4
--- /dev/null
+++ b/doc/theory_manual/Workflows/HERON_Standard_Workflow.md
@@ -0,0 +1,138 @@
+
+## HERON Modes for `standard` Workflow
+For RAVEN-runs-RAVEN (`standard`) mode, users can select:
+- `Opt` Mode
+- `Sweep` Mode
+
+Additionally, users can overwrite either of these modes with:
+- `Debug` Mode
+  - runs limited case and produces diagrams for given IES configuration
+
+
+![HERON_bilevel_optimization](../diagrams/HERON_bilevel_optimization_full.png)
+
+## Outer Level
+In the outer level, regardless of HERON mode (`opt` or `sweep`), we calculate some statistic of a set of economic metrics:
+
+```math
+\mathcal{S}( \mathbb{O} )
+```
+
+where $\mathcal{S}$ is the statistical metric and $\mathbb{O}$ is a set of economic metrics.
+
+### *Sweep Mode*
+
+If selecting sweep mode, the outer level runs a parametric sweep over a list of capacity values $c$ as:
+```math
+\underset{\mathbf{c}}{\text{sweep}} \ \ \mathcal{S}( \mathbb{O} )
+```
+
+Multiple statistics can be reported back per iteration within the parametric sweep, shown here for $K$ different statistics (e.g., expected value, minimum, 95% percentile):
+```math
+\underset{\mathbf{c}}{\text{sweep}} \ \ \Big[\mathcal{S}_0( \mathbb{O} ), \ \mathcal{S}_1( \mathbb{O} ), \ ..., \ \mathcal{S}_K( \mathbb{O} ) \Big]
+```
+
+Additionally, and only in sweep mode, HERON can report back multiple economic metrics, shown below for $M$ different metrics:
+
+```math
+\underset{\mathbf{c}}{\text{sweep}} \ \ \begin{bmatrix} \mathcal{S}_0( \mathbb{O}_0 ), \ \mathcal{S}_1( \mathbb{O}_0 ), \ ..., \ \mathcal{S}_K( \mathbb{O}_0 ) \\ \vdots \\ \mathcal{S}_0( \mathbb{O}_M ), \ \mathcal{S}_1( \mathbb{O}_M ), \ ..., \ \mathcal{S}_K( \mathbb{O}_M ) \end{bmatrix}
+```
+
+NOTE: multiple metrics is not available when selecting levelized cost. Selection for other metrics is limited to NPV, IRR, and PI.
+
+### *Opt Mode*
+
+If selecting opt mode, the outer level uses an optimizer algorithm to determine a set of capacities which optimize (either maximize or minimize) some objective function. Currently, the available optimization algorithms are:
+
+1) Gradient Descent
+2) Bayesian Optimization (default)
+
+Regardless of optimization method, the overall problem is stated as
+
+```math
+\underset{\mathbf{c}}{\text{opt}} \ \ \mathcal{S}( \mathbb{O}_m )
+```
+
+where:
+
+- $\mathbf{c}$ : component (e.g., nuclear, wind turbines) capacities
+- $\mathcal{S}$ : statistical operator (e.g., expected value, value at risk)
+- $\mathbb{O}_m$ : a set of economic metrics from individual inner realizations
+- $m$ : subscript for some economic metric (e.g., net present value, levelized cost)
+- $\mathbf{c}^\star$ : optimal set of capacities for the given objective
+
+Depending on the statistical metric (and even moreso the economic metric), the direction of optimization changes. For example if the economic metric is the net present value, if the statistical metric is expected value then we would want to maximize:
+```math
+\underset{\mathbf{c}}{\text{max}} \ \ \mathbb{E}( \mathbb{O}_{NPV} )
+```
+
+If the statistic is the value-at-risk of the distribution at a given probability level, we'd want to minimize:
+```math
+\underset{\mathbf{c}}{\text{min}} \ \ \mathbb{V}_{0.05}( \mathbb{O}_{NPV} )
+```
+
+The optimal set of capacities for the given objective is
+
+```math
+\mathbf{c}^\star = \underset{\mathbf{c}}{\text{argmax}} \ \mathcal{S}( \mathbb{O}_{NPV} )
+```
+
+## Inner Level
+
+
+We build up that set of economic metrics for a set of possible scenarios/realizations.
+
+$\mathbb{O}_m = \Big[ \ \mathcal{T}({}^s\mathbb{F}) \ :  \ s \in \mathbb{S} \Big]$
+
+- $\mathbb{S}$ : set of scenario indeces, size $S$
+- $s$ : given scenario from $\mathbb{S}$
+- ${}^s\mathbb{F}$ : set of cashflows per scenario (collected over all components, all years, all hours)
+- $\mathcal{T}( )$ : evaluation of a function resulting in an economic metric (typically from a plugin called TEAL)
+
+where $\mathbb{S}$ is a set of scenarios, ${}^s\mathbb{F}$ is a set of cashflows per scenario, and $\mathcal{T}( )$ is an evaluated economic metric (from TEAL). The typical is...
+
+$\mathbb{O}_{NPV} = \Big[ \ NPV({}^s\mathbb{F}) \ : \ s \in \mathbb{S} \Big]$
+
+### for a given scenario ($s$)
+
+the set of all cashflows for that scenario is
+
+${}^s\mathbb{F} = \Bigg[ \bigg[ \Big[{}^s\mathbb{F}_{p,F,i} : p \in \mathbb{P} \Big]: F \in \big[A, B, \Gamma, H \big]\bigg] \ : i \in [+,-] \Bigg ]$
+
+or in other words it is a collection of all cashflow types $F$, either positive or negative, for all components. We can define the set by all cashflow types (each are indexed by years in different ways):
+
+```math
+A_p^i \equiv {}^s\mathbb{F}_{p,i}|_{{}_{F=A}}
+```
+
+```math
+\Gamma_{p,y}^i \equiv {}^s\mathbb{F}_{p,y,i}|_{{}_{F=\Gamma}}
+```
+
+etc.
+
+
+```math
+\mathbb{{}^0F} = \Bigg[ \bigg[ A^i_{p}, \ \Big[ B^i_{p,y}, \ \Gamma^i_{p,y}, \ \sum\limits_{\mathbb{U}_y} \sum\limits_{\mathbb{T}_u} H^i_{p,y}   \ :  \ y \in \mathbb{Y}  \Big] \ :  \ p \in \mathbb{P} \  \bigg ] \ i \in [+,-] \Bigg ]
+```
+
+then the NPV is a linear combination of all cashflows
+
+```math
+NPV(\mathbb{{}^0F} ) = \sum\limits_{{}^0F \in \mathbb{{}^0F}} {}^0F = \sum\limits_{p \in \mathbb{P}}\sum\limits_{y \in \mathbb{Y}} {}^0F_{p,y} = \sum\limits_{p \in \mathbb{P}}\sum\limits_{y \in \mathbb{Y}} {}^0F_{p,y}
+```
+
+and each cashflow is
+
+```math
+{}^0F_{p,y} = \xi_{p,y} \alpha_{p,y} \Big(\frac{d_{p,y}}{d^\prime_{p,y}}\Big)^{\chi_{p,y}}
+```
+
+where
+
+```math
+\xi_{p,y} = \frac{1-\epsilon_{p,F}^\tau \tau}{(1+r)^y (1+\epsilon_{p,F}^\nu \nu)^y} = \frac{1-\epsilon_{p,F}^\tau \tau}{(1+r+\epsilon_{p,F}^\nu \nu + r\epsilon_{p,F}^\nu \nu)^y}
+```
+
+where the $\epsilon$ terms are just a boolean, to toggle if the component cashflow is meant to be taxed ($\tau$) and/or have annualized inflation $\nu$. Each cashflow has an annualized discount rate $r$.
+
diff --git a/doc/theory_manual/diagrams/HERON_bilevel_optimization_full.png b/doc/theory_manual/diagrams/HERON_bilevel_optimization_full.png
new file mode 100644
index 00000000..46d6ad5f
Binary files /dev/null and b/doc/theory_manual/diagrams/HERON_bilevel_optimization_full.png differ
diff --git a/doc/theory_manual/diagrams/HERON_cfs.png b/doc/theory_manual/diagrams/HERON_cfs.png
new file mode 100644
index 00000000..9d7a207f
Binary files /dev/null and b/doc/theory_manual/diagrams/HERON_cfs.png differ
diff --git a/doc/theory_manual/diagrams/HERON_comps.png b/doc/theory_manual/diagrams/HERON_comps.png
new file mode 100644
index 00000000..d01176e4
Binary files /dev/null and b/doc/theory_manual/diagrams/HERON_comps.png differ
diff --git a/doc/theory_manual/diagrams/HERON_time.png b/doc/theory_manual/diagrams/HERON_time.png
new file mode 100644
index 00000000..1c054218
Binary files /dev/null and b/doc/theory_manual/diagrams/HERON_time.png differ