diff --git a/docs/Project.toml b/docs/Project.toml
index a15dc604d..8174017e2 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -19,6 +19,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MarkdownTables = "1862ce21-31c7-451e-824c-f20fa3f90fa2"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Ribasim = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635"
+SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
 
 [compat]
 Configurations = "0.17"
@@ -33,4 +34,5 @@ Legolas = "0.5"
 Logging = "<0.0.1,1"
 MarkdownTables = "1"
 OrderedCollections = "1.6"
+SQLite = "1.5.1"
 julia = "1.10"
diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd
index f80c79a9c..14b16167e 100644
--- a/docs/core/allocation.qmd
+++ b/docs/core/allocation.qmd
@@ -2,15 +2,25 @@
 title: "Allocation"
 ---
 
-Allocation is the process of assigning an allocated abstraction flow rate to user nodes in the model based on information about sources, user demands over various priorities, constraints introduced by nodes, local water availability and graph topology. The allocation procedure implemented in Ribasim is heavily inspired by the [maximum flow problem](https://en.wikipedia.org/wiki/Maximum_flow_problem).
+# Introduction
+Allocation is the process of assigning an allocated abstraction flow rate to user nodes in the physical layer of the model based on information about sources, user demands over various priorities, constraints introduced by nodes, local water availability and graph topology. The allocation procedure implemented in Ribasim is heavily inspired by the [maximum flow problem](https://en.wikipedia.org/wiki/Maximum_flow_problem).
 
-The allocation problem is solved per subnetwork of the Ribasim model. The subnetwork is used to formulate an optimization problem with the [JuMP](https://jump.dev/JuMP.jl/stable/) package, which is solved using the [HiGHS solver](https://highs.dev/). For more information see also the example of solving the maximum flow problem with `JuMP.jl` [here](https://jump.dev/JuMP.jl/stable/tutorials/linear/network_flows/#The-max-flow-problem).
+The allocation problem is solved per subnetwork of the Ribasim model. The subnetwork is used to formulate an optimization problem with the [JuMP](https://jump.dev/JuMP.jl/stable/) package, which is solved using the [HiGHS solver](https://highs.dev/). For more in-depth information see also the example of solving the maximum flow problem with `JuMP.jl` [here](https://jump.dev/JuMP.jl/stable/tutorials/linear/network_flows/#The-max-flow-problem).
 
-# The allocation problem
+# The high level algorithm
+The allocation algorithm consists of 3 types of optimization problems:
+
+1. **Subnetwork demand collection**: Collect the demands of a subnetwork from the main network by optimizing with unlimited capacity from the main network;
+2. **Main network allocation**: Allocate to subnetworks with the above collected demand, and users in the main network;
+3. **Subnetwork allocation**: Allocate to users in the subnetworks with the flows allocated to the subnetwork in the main network allocation.
+
+The total allocation algorithm consists of performing 1 for all subnetworks, then performing 2, then performing 3 for all subnetworks. Not having a main network is also supported, then 1 and 2 are skipped.
+
+# Elements of allocation
 
 The following data of the parameters and state of a Ribasim model are relevant for the allocation problem.
 
-## Allocation problem input
+## Schematisation input
 
 ### The subnetwork
 
@@ -22,7 +32,7 @@ Sources are indicated by a set of edges in the subnetwork
 $$
 E_S^\text{source} \subset \left(S \times S\right) \cap E.
 $$
-That is, if $(i,j) \in E_S^\text{source}$, then $Q_{ij}$ (see the [formal model description](equations.qmd#formal-model-description)) is treated as a source flow in the allocation problem.
+That is, if $(i,j) \in E_S^\text{source}$, then $Q_{ij}$ (see the [formal model description](equations.qmd#formal-model-description)) is treated as a source flow in the allocation problem. These edges are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the main network to a subnetwork.
 
 ### User demands
 
@@ -35,6 +45,8 @@ $$
 On this page we assume that the priorities are given by all integers from $1$ to some $p_{\max} \in \mathbb{N}$. However, in the Ribasim input this is not a requirement; some of these in between priority values can be missing, only the ordering of the given priorities is taken into account.
 :::
 
+## Simulation (physical layer) input
+
 ### Vertical fluxes and local storage
 
 Apart from the source flows denoted by edges, there are other sources of water in the subnetwork, associated with the basins in the subnetwork $B_S = B \cap S$. Firstly there is the sum of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each basin
@@ -48,20 +60,20 @@ $$
 $$
 Note that this value can be negative, which we interpret as a demand from the basin.
 
-### Flow magnitude and direction constraints
+### Constraining factors
 
+#### Flow magnitude and direction constraints
 Nodes in the Ribasim model that have a `max_flow_rate`, i.e. pumps and outlets, put a constraint on the flow through that node. Some nodes only allow flow in one direction, like pumps, outlets and tabulated rating curves.
 
-### Fractional flows and user return flows
-
+#### Fractional flows and user return flows
 Both fractional flow nodes and user nodes dictate proportional relationships between flows over edges in the subnetwork. Users have a return factor $0 \le r_i \le 1, i \in U_S$.
 
-## The allocation optimization problem
-
-### The allocation subgraph
+## The allocation graph
 
 A new graph is created from the subnetwork, which we call an allocation graph. The allocation graph is almost a subgraph of the main (flow) model, apart from the fact that an allocation graph can contain edges which are not in the main model.
 
+### Nodes and edges
+
 The allocation graph consists of:
 
 - Nodes $V'_S \subset V_S$, where each basin, source and user in the subnetwork get a node in the allocation graph. Also nodes that have fractional flow outneighbors get a node in the allocation graph.
@@ -74,11 +86,11 @@ For notational convenience, we use the notation
     V^{\text{in}}_S(j) = \left\{i \in V'_S : (i,j) \in E_S\right\}
 \end{align}
 
-for the set of in-neighbors and out-neighbors of a node in the allocation graph respectively
+for the set of in-neighbors and out-neighbors of a node in the allocation graph respectively.
 
-### The allocation graph capacities
+### Capacities
 
-The capacities of the edges of the allocation graph are collected in the sparse capacity matrix $C_S \in \overline{\mathbb{R}}_{\ge 0}^{n'\times n'}$ where $n' = \#V'_S$ is the number of nodes in the allocation graph. The capacities can be infinite.
+Each edge in the allocation graph has an associated capacity. These capacities are collected in the sparse capacity matrix $C_S \in \overline{\mathbb{R}}_{\ge 0}^{n'\times n'}$ where $n' = \#V'_S$ is the number of nodes in the allocation graph. The capacities can be infinite, if there is nothing in the model constraining the capacity of the edge.
 
 The capacities are determined in 4 different ways:
 
@@ -86,9 +98,13 @@ The capacities are determined in 4 different ways:
 - The capacity of the edge $e \in E_S$ is given by the smallest `max_flow_rate` of the nodes along the equivalent edges in the subnetwork. If there are no nodes with a `max_flow_rate`, the edge capacity is infinite;
 - If the edge is a source, the capacity of the edge is given by the flow rate of that source.
 
-### The optimization variables
+# The optimization problem
 
-There are 2 types of variables whose value has to be determined to solve the allocation problem:
+The optimization problem for a subnetwork is a linear optimization problem consisting of an objective function with associated constraints on a set of variables, all of which are introduced below.
+
+## The optimization variables
+
+There are 2 types of variable whose value has to be determined to solve the allocation problem:
 
 - The flows $F \in \mathbb{R}_{\ge 0}^{n'\times n'}$ over the edges in the allocation graph;
 - The allocations to the basins
@@ -100,7 +116,7 @@ $$
 Currently the basin allocations are not taken into account in the implementation.
 :::
 
-### The optimization objective
+## The optimization objective
 
 The goal of allocation is to get the flow to the users as close as possible to their demand. To achieve this, the following objectives are supported:
 
@@ -121,6 +137,11 @@ $$
     \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right| + c \sum_{e \in E_S} F_e
 $$
 
+:::{.callout-note}
+When performing main network allocation, the connections to subnetworks are also interpreted as users with demands determined by subnetwork demand collection.
+:::
+
+
 To avoid division by $0$ errors, if a `*_relative` objective is used and a demand is $0$, the coefficient of the flow $F_{ij}$ is set to $0$.
 
 For `*_absolute` objectives the optimizer cares about the actual amount of water allocated to a user, for `*_relative` objectives it cares about the fraction of the demand allocated to the user. For `quadratic_*` objectives the optimizer cares about avoiding large shortages, for `linear_*` objectives it treats all deviations equally.
@@ -137,7 +158,7 @@ The absolute value applied here is not supported in a linear programming context
 In the future new optimization objectives will be introduced, for demands of basins and priorities over sources. These will be used in combination with the above, in the form of goal programming.
 :::
 
-### The optimization variable constraints
+## The optimization constraints
 - Flow conservation: For the basins in the allocation graph we have that
 $$
     \sum_{j=1}^{n'} F_{kj} \le \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S.
@@ -148,6 +169,12 @@ $$
     F_{ij} \le \left(C_S\right)_{ij}, \quad \forall(i,j) \in E_S.
 $$ {#eq-capacityconstraint}
 By the definition of $C_S$ this also includes the source flows.
+
+:::{.callout-note}
+When performing subnetwork demand collection, these capacities are set to $\infty$ for edges which connect the main network to a subnetwork.
+:::
+
+
 - User outflow: The outflow of the user is dictated by the inflow and the return factor:
 $$
     F_{ik} = r_k \cdot F_{kj} \quad
@@ -200,6 +227,26 @@ In the future there will be 2 more optimization solves:
 
 ## Example
 
-:::{.callout-note}
-An example with figures and data will be added here after addition of [allocation output files](https://github.com/Deltares/Ribasim/issues/659).
-:::
+The following is an example of an optimization problem for the example shown [here](../python/examples.ipynb#model-with-allocation):
+
+
+```{julia}
+# | code-fold: true
+using Ribasim
+using SQLite
+
+toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml")
+p = Ribasim.Model(toml_path).integrator.p
+
+allocation_model = p.allocation_models[1]
+t = 0.0
+priority_idx = 1
+
+Ribasim.set_flow!(p.graph, Ribasim.NodeID(1), Ribasim.NodeID(2), 1.0)
+
+Ribasim.adjust_source_flows!(allocation_model, p, priority_idx)
+Ribasim.adjust_edge_capacities!(allocation_model, p, priority_idx)
+Ribasim.set_objective_priority!(allocation_model, p, t, priority_idx)
+
+println(p.allocation_models[1].problem)
+```
diff --git a/pixi.toml b/pixi.toml
index 1f1d5183e..8732b2483 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -41,11 +41,11 @@ build-julia-docs = { cmd = "julia --project=docs docs/make.jl", depends_on = [
 ] }
 quartodoc-build = { cmd = "quartodoc build && rm objects.json", cwd = "docs" }
 quarto-preview = { cmd = "quarto preview docs", depends_on = [
-    "quartodoc-build",
+    "quartodoc-build", "generate-testmodels"
 ] }
 quarto-check = { cmd = "quarto check all", depends_on = ["quartodoc-build"] }
 quarto-render = { cmd = "julia --project=docs --eval 'using Pkg; Pkg.build(\"IJulia\")' && quarto render docs --to html --execute", depends_on = [
-    "quartodoc-build",
+    "quartodoc-build", "generate-testmodels"
 ] }
 docs = { depends_on = ["build-julia-docs", "quarto-preview"] }
 # Lint