diff --git a/docs/contracts/known-problems.md b/docs/contracts/known-problems.md
index a29ac9c1..0dfc30e1 100644
--- a/docs/contracts/known-problems.md
+++ b/docs/contracts/known-problems.md
@@ -15,23 +15,39 @@ mathematical formulation, a runnable Arco model, and expected solution values.
 
 Known optimum: `x = 1`, `y = 4`, objective `11`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import arco
+model = arco.Model()
+x = model.add_variable(bounds=arco.Bounds(lower=1.0, upper=float("inf")), name="x")
+y = model.add_variable(bounds=arco.Bounds(lower=2.0, upper=float("inf")), name="y")
+model.add_constraint(x + y >= 5.0, name="demand")
+model.minimize(3.0 * x + 2.0 * y)
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+round(solution.objective_value, 6)
+round(solution.get_primal(index=x), 6)
+round(solution.get_primal(index=y), 6)
+```
+
 ```python doctest
->>> import arco
->>> model = arco.Model()
->>> x = model.add_variable(bounds=arco.Bounds(lower=1.0, upper=float("inf")), name="x")
->>> y = model.add_variable(bounds=arco.Bounds(lower=2.0, upper=float("inf")), name="y")
->>> model.add_constraint(x + y >= 5.0, name="demand")
-Constraint('demand', Bounds(5, inf))
->>> model.minimize(3.0 * x + 2.0 * y)
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
-11.0
->>> round(solution.get_primal(index=x), 6)
-1.0
->>> round(solution.get_primal(index=y), 6)
-4.0
+# >>> import arco
+# >>> model = arco.Model()
+# >>> x = model.add_variable(bounds=arco.Bounds(lower=1.0, upper=float("inf")), name="x")
+# >>> y = model.add_variable(bounds=arco.Bounds(lower=2.0, upper=float("inf")), name="y")
+# >>> model.add_constraint(x + y >= 5.0, name="demand")
+# Constraint('demand', Bounds(5, inf))
+# >>> model.minimize(3.0 * x + 2.0 * y)
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# 11.0
+# >>> round(solution.get_primal(index=x), 6)
+# 1.0
+# >>> round(solution.get_primal(index=y), 6)
+# 4.0
 ```
 
 ## Production Mix
@@ -47,27 +63,45 @@ True
 
 Known optimum: `x = 20`, `y = 60`, objective `2600`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import arco
+model = arco.Model()
+x = model.add_variable(bounds=arco.NonNegativeFloat, name="x")
+y = model.add_variable(bounds=arco.NonNegativeFloat, name="y")
+model.add_constraint(x <= 40.0, name="demand")
+model.add_constraint(x + y <= 80.0, name="labor_a")
+model.add_constraint(2.0 * x + y <= 100.0, name="labor_b")
+model.maximize(40.0 * x + 30.0 * y)
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+round(solution.objective_value, 6)
+round(solution.get_value(x), 6)
+round(solution.get_value(y), 6)
+```
+
 ```python doctest
->>> import arco
->>> model = arco.Model()
->>> x = model.add_variable(bounds=arco.NonNegativeFloat, name="x")
->>> y = model.add_variable(bounds=arco.NonNegativeFloat, name="y")
->>> model.add_constraint(x <= 40.0, name="demand")
-Constraint('demand', Bounds(-inf, 40))
->>> model.add_constraint(x + y <= 80.0, name="labor_a")
-Constraint('labor_a', Bounds(-inf, 80))
->>> model.add_constraint(2.0 * x + y <= 100.0, name="labor_b")
-Constraint('labor_b', Bounds(-inf, 100))
->>> model.maximize(40.0 * x + 30.0 * y)
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
-2600.0
->>> round(solution.get_value(x), 6)
-20.0
->>> round(solution.get_value(y), 6)
-60.0
+# >>> import arco
+# >>> model = arco.Model()
+# >>> x = model.add_variable(bounds=arco.NonNegativeFloat, name="x")
+# >>> y = model.add_variable(bounds=arco.NonNegativeFloat, name="y")
+# >>> model.add_constraint(x <= 40.0, name="demand")
+# Constraint('demand', Bounds(-inf, 40))
+# >>> model.add_constraint(x + y <= 80.0, name="labor_a")
+# Constraint('labor_a', Bounds(-inf, 80))
+# >>> model.add_constraint(2.0 * x + y <= 100.0, name="labor_b")
+# Constraint('labor_b', Bounds(-inf, 100))
+# >>> model.maximize(40.0 * x + 30.0 * y)
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# 2600.0
+# >>> round(solution.get_value(x), 6)
+# 20.0
+# >>> round(solution.get_value(y), 6)
+# 60.0
 ```
 
 ## $A x = b$ Feasibility
@@ -87,27 +121,47 @@ A=\begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix},\quad b=\begin{bmatrix}3 \\ 1\end{b
 
 Known solution: `x = [1, 2]`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import numpy as np
+import arco
+A = np.array([[1.0, 1.0], [1.0, 0.0]], dtype=float)
+b = np.array([3.0, 1.0], dtype=float)
+model = arco.Model()
+j = arco.IndexSet(name="j", size=2)
+x = model.add_variables(j, bounds=arco.NonNegativeFloat, name="x")
+row_constraints = model.add_constraints((A @ x) == b, name="row")
+len(row_constraints)
+row_constraints[0].bounds
+row_constraints[1].bounds
+model.minimize(0.0)
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+tuple(float(v) for v in np.round(solution.get_value(x), 6))
+```
+
 ```python doctest
->>> import numpy as np
->>> import arco
->>> A = np.array([[1.0, 1.0], [1.0, 0.0]], dtype=float)
->>> b = np.array([3.0, 1.0], dtype=float)
->>> model = arco.Model()
->>> j = arco.IndexSet(name="j", size=2)
->>> x = model.add_variables(j, bounds=arco.NonNegativeFloat, name="x")
->>> row_constraints = model.add_constraints((A @ x) == b, name="row")
->>> len(row_constraints)
-2
->>> row_constraints[0].bounds
-Bounds(lower=3, upper=3)
->>> row_constraints[1].bounds
-Bounds(lower=1, upper=1)
->>> model.minimize(0.0)
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> tuple(float(v) for v in np.round(solution.get_value(x), 6))
-(1.0, 2.0)
+# >>> import numpy as np
+# >>> import arco
+# >>> A = np.array([[1.0, 1.0], [1.0, 0.0]], dtype=float)
+# >>> b = np.array([3.0, 1.0], dtype=float)
+# >>> model = arco.Model()
+# >>> j = arco.IndexSet(name="j", size=2)
+# >>> x = model.add_variables(j, bounds=arco.NonNegativeFloat, name="x")
+# >>> row_constraints = model.add_constraints((A @ x) == b, name="row")
+# >>> len(row_constraints)
+# 2
+# >>> row_constraints[0].bounds
+# Bounds(lower=3, upper=3)
+# >>> row_constraints[1].bounds
+# Bounds(lower=1, upper=1)
+# >>> model.minimize(0.0)
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> tuple(float(v) for v in np.round(solution.get_value(x), 6))
+# (1.0, 2.0)
 ```
 
 ## 0-1 Knapsack
@@ -123,26 +177,46 @@ True
 With values `[6,5,4]` and weights `[3,2,1]`, the known optimum value is `10`
 (select items 1 and 3).
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import arco
+items = ["a", "b", "c"]
+values = [6.0, 5.0, 4.0]
+weights = [3.0, 2.0, 1.0]
+model = arco.Model()
+x = {name: model.add_variable(bounds=arco.Binary, name=name) for name in items}
+model.add_constraint(
+    sum(w * x[name] for w, name in zip(weights, items, strict=True)) <= 4.0,
+    name="capacity",
+)
+model.maximize(sum(v * x[name] for v, name in zip(values, items, strict=True)))
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+round(solution.objective_value, 6)
+tuple(round(solution.get_value(x[name]), 6) for name in items)
+```
+
 ```python doctest
->>> import arco
->>> items = ["a", "b", "c"]
->>> values = [6.0, 5.0, 4.0]
->>> weights = [3.0, 2.0, 1.0]
->>> model = arco.Model()
->>> x = {name: model.add_variable(bounds=arco.Binary, name=name) for name in items}
->>> model.add_constraint(
-...     sum(w * x[name] for w, name in zip(weights, items, strict=True)) <= 4.0,
-...     name="capacity",
-... )
-Constraint('capacity', Bounds(-inf, 4))
->>> model.maximize(sum(v * x[name] for v, name in zip(values, items, strict=True)))
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
-10.0
->>> tuple(round(solution.get_value(x[name]), 6) for name in items)
-(1.0, 0.0, 1.0)
+# >>> import arco
+# >>> items = ["a", "b", "c"]
+# >>> values = [6.0, 5.0, 4.0]
+# >>> weights = [3.0, 2.0, 1.0]
+# >>> model = arco.Model()
+# >>> x = {name: model.add_variable(bounds=arco.Binary, name=name) for name in items}
+# >>> model.add_constraint(
+# ...     sum(w * x[name] for w, name in zip(weights, items, strict=True)) <= 4.0,
+# ...     name="capacity",
+# ... )
+# Constraint('capacity', Bounds(-inf, 4))
+# >>> model.maximize(sum(v * x[name] for v, name in zip(values, items, strict=True)))
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# 10.0
+# >>> tuple(round(solution.get_value(x[name]), 6) for name in items)
+# (1.0, 0.0, 1.0)
 ```
 
 ## Network Flow (numpy integration)
@@ -159,31 +233,54 @@ True
 Known optimum on the small test graph: objective `2` using the direct arc
 `0 \to 1`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import numpy as np
+import arco
+costs = np.array([[0.0, 2.0, 5.0], [0.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+balance = np.array([1.0, -1.0, 0.0])
+model = arco.Model()
+src = arco.IndexSet(name="src", size=3)
+dst = arco.IndexSet(name="dst", size=3)
+x = model.add_variables(src, dst, bounds=arco.Binary, name="x")
+no_arc = model.add_constraints(x[costs == 0.0] == 0.0, name="no_arc")
+len(no_arc)
+flow = model.add_constraints((x.sum(over=dst) - x.sum(over=src)) == balance, name="flow")
+len(flow)
+model.minimize(np.sum(costs * x))
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+round(solution.objective_value, 6)
+round(solution.get_value(x[0, 1]), 6)
+abs(solution.get_value(x[0, 2])) < 1e-6
+```
+
 ```python doctest
->>> import numpy as np
->>> import arco
->>> costs = np.array([[0.0, 2.0, 5.0], [0.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
->>> balance = np.array([1.0, -1.0, 0.0])
->>> model = arco.Model()
->>> src = arco.IndexSet(name="src", size=3)
->>> dst = arco.IndexSet(name="dst", size=3)
->>> x = model.add_variables(src, dst, bounds=arco.Binary, name="x")
->>> no_arc = model.add_constraints(x[costs == 0.0] == 0.0, name="no_arc")
->>> len(no_arc)
-6
->>> flow = model.add_constraints((x.sum(over=dst) - x.sum(over=src)) == balance, name="flow")
->>> len(flow)
-3
->>> model.minimize(np.sum(costs * x))
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
-2.0
->>> round(solution.get_value(x[0, 1]), 6)
-1.0
->>> abs(solution.get_value(x[0, 2])) < 1e-6
-True
+# >>> import numpy as np
+# >>> import arco
+# >>> costs = np.array([[0.0, 2.0, 5.0], [0.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+# >>> balance = np.array([1.0, -1.0, 0.0])
+# >>> model = arco.Model()
+# >>> src = arco.IndexSet(name="src", size=3)
+# >>> dst = arco.IndexSet(name="dst", size=3)
+# >>> x = model.add_variables(src, dst, bounds=arco.Binary, name="x")
+# >>> no_arc = model.add_constraints(x[costs == 0.0] == 0.0, name="no_arc")
+# >>> len(no_arc)
+# 6
+# >>> flow = model.add_constraints((x.sum(over=dst) - x.sum(over=src)) == balance, name="flow")
+# >>> len(flow)
+# 3
+# >>> model.minimize(np.sum(costs * x))
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# 2.0
+# >>> round(solution.get_value(x[0, 1]), 6)
+# 1.0
+# >>> abs(solution.get_value(x[0, 2])) < 1e-6
+# True
 ```
 
 ## Unit Commitment
@@ -200,32 +297,56 @@ True
 Known optimum for this two-unit instance: objective `275`, generation
 `[100, 20]`, commitment `[1, 1]`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import numpy as np
+import arco
+gen_max = np.array([100.0, 80.0])
+fixed_cost = np.array([10.0, 5.0])
+variable_cost = np.array([2.0, 3.0])
+demand = 120.0
+model = arco.Model()
+units = arco.IndexSet(name="unit", size=2)
+g = model.add_variables(units, bounds=arco.NonNegativeFloat, name="gen")
+u = model.add_variables(units, bounds=arco.Binary, name="commit")
+model.add_constraint(g.sum() == demand, name="balance")
+capacity = model.add_constraints(g <= gen_max * u, name="capacity")
+len(capacity)
+model.minimize(np.sum(fixed_cost * u + variable_cost * g))
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+round(solution.objective_value, 6)
+tuple(float(v) for v in np.round(solution.get_value(g), 6))
+tuple(float(v) for v in np.round(solution.get_value(u), 6))
+```
+
 ```python doctest
->>> import numpy as np
->>> import arco
->>> gen_max = np.array([100.0, 80.0])
->>> fixed_cost = np.array([10.0, 5.0])
->>> variable_cost = np.array([2.0, 3.0])
->>> demand = 120.0
->>> model = arco.Model()
->>> units = arco.IndexSet(name="unit", size=2)
->>> g = model.add_variables(units, bounds=arco.NonNegativeFloat, name="gen")
->>> u = model.add_variables(units, bounds=arco.Binary, name="commit")
->>> model.add_constraint(g.sum() == demand, name="balance")
-Constraint('balance', Bounds(120, 120))
->>> capacity = model.add_constraints(g <= gen_max * u, name="capacity")
->>> len(capacity)
-2
->>> model.minimize(np.sum(fixed_cost * u + variable_cost * g))
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
-275.0
->>> tuple(float(v) for v in np.round(solution.get_value(g), 6))
-(100.0, 20.0)
->>> tuple(float(v) for v in np.round(solution.get_value(u), 6))
-(1.0, 1.0)
+# >>> import numpy as np
+# >>> import arco
+# >>> gen_max = np.array([100.0, 80.0])
+# >>> fixed_cost = np.array([10.0, 5.0])
+# >>> variable_cost = np.array([2.0, 3.0])
+# >>> demand = 120.0
+# >>> model = arco.Model()
+# >>> units = arco.IndexSet(name="unit", size=2)
+# >>> g = model.add_variables(units, bounds=arco.NonNegativeFloat, name="gen")
+# >>> u = model.add_variables(units, bounds=arco.Binary, name="commit")
+# >>> model.add_constraint(g.sum() == demand, name="balance")
+# Constraint('balance', Bounds(120, 120))
+# >>> capacity = model.add_constraints(g <= gen_max * u, name="capacity")
+# >>> len(capacity)
+# 2
+# >>> model.minimize(np.sum(fixed_cost * u + variable_cost * g))
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# 275.0
+# >>> tuple(float(v) for v in np.round(solution.get_value(g), 6))
+# (100.0, 20.0)
+# >>> tuple(float(v) for v in np.round(solution.get_value(u), 6))
+# (1.0, 1.0)
 ```
 
 ## Network Design
@@ -243,43 +364,76 @@ True
 Known optimum for this toy graph: select edges `0->1` and `1->3`, objective
 `-4.8`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import numpy as np
+import arco
+G = np.array([
+    [0.0, 5.0, 4.0, 0.0],
+    [0.0, 0.0, 0.0, 5.0],
+    [0.0, 0.0, 0.0, 4.0],
+    [0.0, 0.0, 0.0, 0.0],
+])
+model = arco.Model()
+src = arco.IndexSet(name="src", size=4)
+dst = arco.IndexSet(name="dst", size=4)
+edge = model.add_variables(src, dst, bounds=arco.Binary, name="edge")
+flow = model.add_variables(src, dst, bounds=arco.NonNegativeFloat, name="flow")
+model.add_constraint(np.sum(edge) <= 2.0, name="edge_budget")
+capacity = model.add_constraints(flow <= G * edge, name="capacity")
+len(capacity)
+conservation = model.add_constraints(
+    flow[1:-1, :].sum(over=dst) == flow[:, 1:-1].sum(over=src),
+    name="flow_conservation",
+)
+len(conservation)
+model.minimize(0.1 * np.sum(edge) - flow[0, :].sum())
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+round(solution.objective_value, 6)
+solution.get_value(edge[0, 1]) > 0.9
+solution.get_value(edge[1, 3]) > 0.9
+abs(solution.get_value(edge[0, 2])) < 1e-6
+```
+
 ```python doctest
->>> import numpy as np
->>> import arco
->>> G = np.array([
-...     [0.0, 5.0, 4.0, 0.0],
-...     [0.0, 0.0, 0.0, 5.0],
-...     [0.0, 0.0, 0.0, 4.0],
-...     [0.0, 0.0, 0.0, 0.0],
-... ])
->>> model = arco.Model()
->>> src = arco.IndexSet(name="src", size=4)
->>> dst = arco.IndexSet(name="dst", size=4)
->>> edge = model.add_variables(src, dst, bounds=arco.Binary, name="edge")
->>> flow = model.add_variables(src, dst, bounds=arco.NonNegativeFloat, name="flow")
->>> model.add_constraint(np.sum(edge) <= 2.0, name="edge_budget")
-Constraint('edge_budget', Bounds(-inf, 2))
->>> capacity = model.add_constraints(flow <= G * edge, name="capacity")
->>> len(capacity)
-16
->>> conservation = model.add_constraints(
-...     flow[1:-1, :].sum(over=dst) == flow[:, 1:-1].sum(over=src),
-...     name="flow_conservation",
-... )
->>> len(conservation)
-2
->>> model.minimize(0.1 * np.sum(edge) - flow[0, :].sum())
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
--4.8
->>> solution.get_value(edge[0, 1]) > 0.9
-True
->>> solution.get_value(edge[1, 3]) > 0.9
-True
->>> abs(solution.get_value(edge[0, 2])) < 1e-6
-True
+# >>> import numpy as np
+# >>> import arco
+# >>> G = np.array([
+# ...     [0.0, 5.0, 4.0, 0.0],
+# ...     [0.0, 0.0, 0.0, 5.0],
+# ...     [0.0, 0.0, 0.0, 4.0],
+# ...     [0.0, 0.0, 0.0, 0.0],
+# ... ])
+# >>> model = arco.Model()
+# >>> src = arco.IndexSet(name="src", size=4)
+# >>> dst = arco.IndexSet(name="dst", size=4)
+# >>> edge = model.add_variables(src, dst, bounds=arco.Binary, name="edge")
+# >>> flow = model.add_variables(src, dst, bounds=arco.NonNegativeFloat, name="flow")
+# >>> model.add_constraint(np.sum(edge) <= 2.0, name="edge_budget")
+# Constraint('edge_budget', Bounds(-inf, 2))
+# >>> capacity = model.add_constraints(flow <= G * edge, name="capacity")
+# >>> len(capacity)
+# 16
+# >>> conservation = model.add_constraints(
+# ...     flow[1:-1, :].sum(over=dst) == flow[:, 1:-1].sum(over=src),
+# ...     name="flow_conservation",
+# ... )
+# >>> len(conservation)
+# 2
+# >>> model.minimize(0.1 * np.sum(edge) - flow[0, :].sum())
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# -4.8
+# >>> solution.get_value(edge[0, 1]) > 0.9
+# True
+# >>> solution.get_value(edge[1, 3]) > 0.9
+# True
+# >>> abs(solution.get_value(edge[0, 2])) < 1e-6
+# True
 ```
 
 ## Energy Storage Dispatch
@@ -297,40 +451,71 @@ True
 Known optimum for `T=2` with initial storage `2`: objective `40`, thermal
 generation `[2,2]`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import numpy as np
+import arco
+T = 2
+demand = np.array([5.0, 5.0])
+available = np.array([2.0, 2.0])
+initial_storage = 2.0
+model = arco.Model()
+time = arco.IndexSet(name="t", size=T)
+r = model.add_variables(time, bounds=arco.NonNegativeFloat, name="renewable")
+p = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=10.0), name="thermal")
+s = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=2.0), name="storage")
+c = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=1.0), name="charge")
+d = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=1.0), name="discharge")
+balance = model.add_constraints(p + r + d == demand + c, name="balance")
+len(balance)
+s_prev = np.concatenate([[initial_storage], s[:-1]])
+dynamics = model.add_constraints(s == s_prev + c - d, name="storage_dynamics")
+len(dynamics)
+renewable_cap = model.add_constraints(r <= available, name="renewable_cap")
+len(renewable_cap)
+model.minimize(10.0 * p.sum())
+solution = model.solve(log_to_console=False)
+solution.is_optimal()
+round(solution.objective_value, 6)
+tuple(float(v) for v in np.round(solution.get_value(p), 6))
+tuple(float(v) for v in np.round(solution.get_value(d), 6))
+```
+
 ```python doctest
->>> import numpy as np
->>> import arco
->>> T = 2
->>> demand = np.array([5.0, 5.0])
->>> available = np.array([2.0, 2.0])
->>> initial_storage = 2.0
->>> model = arco.Model()
->>> time = arco.IndexSet(name="t", size=T)
->>> r = model.add_variables(time, bounds=arco.NonNegativeFloat, name="renewable")
->>> p = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=10.0), name="thermal")
->>> s = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=2.0), name="storage")
->>> c = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=1.0), name="charge")
->>> d = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=1.0), name="discharge")
->>> balance = model.add_constraints(p + r + d == demand + c, name="balance")
->>> len(balance)
-2
->>> s_prev = np.concatenate([[initial_storage], s[:-1]])
->>> dynamics = model.add_constraints(s == s_prev + c - d, name="storage_dynamics")
->>> len(dynamics)
-2
->>> renewable_cap = model.add_constraints(r <= available, name="renewable_cap")
->>> len(renewable_cap)
-2
->>> model.minimize(10.0 * p.sum())
->>> solution = model.solve(log_to_console=False)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
-40.0
->>> tuple(float(v) for v in np.round(solution.get_value(p), 6))
-(2.0, 2.0)
->>> tuple(float(v) for v in np.round(solution.get_value(d), 6))
-(1.0, 1.0)
+# >>> import numpy as np
+# >>> import arco
+# >>> T = 2
+# >>> demand = np.array([5.0, 5.0])
+# >>> available = np.array([2.0, 2.0])
+# >>> initial_storage = 2.0
+# >>> model = arco.Model()
+# >>> time = arco.IndexSet(name="t", size=T)
+# >>> r = model.add_variables(time, bounds=arco.NonNegativeFloat, name="renewable")
+# >>> p = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=10.0), name="thermal")
+# >>> s = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=2.0), name="storage")
+# >>> c = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=1.0), name="charge")
+# >>> d = model.add_variables(time, bounds=arco.Bounds(lower=0.0, upper=1.0), name="discharge")
+# >>> balance = model.add_constraints(p + r + d == demand + c, name="balance")
+# >>> len(balance)
+# 2
+# >>> s_prev = np.concatenate([[initial_storage], s[:-1]])
+# >>> dynamics = model.add_constraints(s == s_prev + c - d, name="storage_dynamics")
+# >>> len(dynamics)
+# 2
+# >>> renewable_cap = model.add_constraints(r <= available, name="renewable_cap")
+# >>> len(renewable_cap)
+# 2
+# >>> model.minimize(10.0 * p.sum())
+# >>> solution = model.solve(log_to_console=False)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# 40.0
+# >>> tuple(float(v) for v in np.round(solution.get_value(p), 6))
+# (2.0, 2.0)
+# >>> tuple(float(v) for v in np.round(solution.get_value(d), 6))
+# (1.0, 1.0)
 ```
 
 ## N-Queens
@@ -348,39 +533,68 @@ True
 
 For `N=4`, the model is feasible with objective `0`.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import numpy as np
+import arco
+N = 4
+model = arco.Model()
+rows = arco.IndexSet(name="row", size=N)
+cols = arco.IndexSet(name="col", size=N)
+q = model.add_variables(rows, cols, bounds=arco.Binary, name="queen")
+row_cons = model.add_constraints(q.sum(over=cols) == 1.0, name="row")
+len(row_cons)
+col_cons = model.add_constraints(q.sum(over=rows) == 1.0, name="col")
+len(col_cons)
+diag_down = [model.add_constraint(np.diag(q, k).sum() <= 1.0, name=f"diag_down_{k}") for k in range(-(N - 1), N)]
+len(diag_down)
+diag_up = [model.add_constraint(np.diag(np.fliplr(q), k).sum() <= 1.0, name=f"diag_up_{k}") for k in range(-(N - 1), N)]
+len(diag_up)
+model.minimize(0.0)
+solution = model.solve(log_to_console=False)
+board = solution.get_value(q).reshape(N, N)
+solution.is_optimal()
+round(solution.objective_value, 6)
+tuple(float(v) for v in np.round(np.sum(board, axis=0), 6))
+tuple(float(v) for v in np.round(np.sum(board, axis=1), 6))
+bool(max(np.diag(board, k).sum() for k in range(-(N - 1), N)) <= 1.0 + 1e-6)
+bool(max(np.diag(np.fliplr(board), k).sum() for k in range(-(N - 1), N)) <= 1.0 + 1e-6)
+```
+
 ```python doctest
->>> import numpy as np
->>> import arco
->>> N = 4
->>> model = arco.Model()
->>> rows = arco.IndexSet(name="row", size=N)
->>> cols = arco.IndexSet(name="col", size=N)
->>> q = model.add_variables(rows, cols, bounds=arco.Binary, name="queen")
->>> row_cons = model.add_constraints(q.sum(over=cols) == 1.0, name="row")
->>> len(row_cons)
-4
->>> col_cons = model.add_constraints(q.sum(over=rows) == 1.0, name="col")
->>> len(col_cons)
-4
->>> diag_down = [model.add_constraint(np.diag(q, k).sum() <= 1.0, name=f"diag_down_{k}") for k in range(-(N - 1), N)]
->>> len(diag_down)
-7
->>> diag_up = [model.add_constraint(np.diag(np.fliplr(q), k).sum() <= 1.0, name=f"diag_up_{k}") for k in range(-(N - 1), N)]
->>> len(diag_up)
-7
->>> model.minimize(0.0)
->>> solution = model.solve(log_to_console=False)
->>> board = solution.get_value(q).reshape(N, N)
->>> solution.is_optimal()
-True
->>> round(solution.objective_value, 6)
-0.0
->>> tuple(float(v) for v in np.round(np.sum(board, axis=0), 6))
-(1.0, 1.0, 1.0, 1.0)
->>> tuple(float(v) for v in np.round(np.sum(board, axis=1), 6))
-(1.0, 1.0, 1.0, 1.0)
->>> bool(max(np.diag(board, k).sum() for k in range(-(N - 1), N)) <= 1.0 + 1e-6)
-True
->>> bool(max(np.diag(np.fliplr(board), k).sum() for k in range(-(N - 1), N)) <= 1.0 + 1e-6)
-True
+# >>> import numpy as np
+# >>> import arco
+# >>> N = 4
+# >>> model = arco.Model()
+# >>> rows = arco.IndexSet(name="row", size=N)
+# >>> cols = arco.IndexSet(name="col", size=N)
+# >>> q = model.add_variables(rows, cols, bounds=arco.Binary, name="queen")
+# >>> row_cons = model.add_constraints(q.sum(over=cols) == 1.0, name="row")
+# >>> len(row_cons)
+# 4
+# >>> col_cons = model.add_constraints(q.sum(over=rows) == 1.0, name="col")
+# >>> len(col_cons)
+# 4
+# >>> diag_down = [model.add_constraint(np.diag(q, k).sum() <= 1.0, name=f"diag_down_{k}") for k in range(-(N - 1), N)]
+# >>> len(diag_down)
+# 7
+# >>> diag_up = [model.add_constraint(np.diag(np.fliplr(q), k).sum() <= 1.0, name=f"diag_up_{k}") for k in range(-(N - 1), N)]
+# >>> len(diag_up)
+# 7
+# >>> model.minimize(0.0)
+# >>> solution = model.solve(log_to_console=False)
+# >>> board = solution.get_value(q).reshape(N, N)
+# >>> solution.is_optimal()
+# True
+# >>> round(solution.objective_value, 6)
+# 0.0
+# >>> tuple(float(v) for v in np.round(np.sum(board, axis=0), 6))
+# (1.0, 1.0, 1.0, 1.0)
+# >>> tuple(float(v) for v in np.round(np.sum(board, axis=1), 6))
+# (1.0, 1.0, 1.0, 1.0)
+# >>> bool(max(np.diag(board, k).sum() for k in range(-(N - 1), N)) <= 1.0 + 1e-6)
+# True
+# >>> bool(max(np.diag(np.fliplr(board), k).sum() for k in range(-(N - 1), N)) <= 1.0 + 1e-6)
+# True
 ```
diff --git a/docs/contracts/python-model-rendering.md b/docs/contracts/python-model-rendering.md
index 08b23d1d..c6675db3 100644
--- a/docs/contracts/python-model-rendering.md
+++ b/docs/contracts/python-model-rendering.md
@@ -2,90 +2,166 @@
 
 This file encodes rendering guarantees for the Python API as executable doctests.
 
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import arco
+model = arco.Model()
+x1 = model.add_variable(bounds=arco.Binary, name="x[1]")
+x2 = model.add_variable(bounds=arco.Binary, name="x[2]")
+x3 = model.add_variable(bounds=arco.Binary, name="x[3]")
+x4 = model.add_variable(bounds=arco.Binary, name="x[4]")
+x5 = model.add_variable(bounds=arco.Binary, name="x[5]")
+_ = model.add_constraint(2.0 * x1 + 8.0 * x2 + 4.0 * x3 + 2.0 * x4 + 5.0 * x5 <= 10.0)
+_ = model.add_constraint(x1 + x3 - x5 == 0.0)
+model.maximize(5.0 * x1 + 3.0 * x2 + 2.0 * x3 + 7.0 * x4 + 4.0 * x5)
+rendered = str(model)
+rendered.startswith("Max 5 x[1] + 3 x[2] + 2 x[3] + 7 x[4] + 4 x[5]")
+"\ns.t.\n" in rendered
+"Subject to" in rendered
+"Binary: x[1], x[2], x[3], x[4], x[5]" in rendered
+" <= 10" in rendered
+operator_columns = [
+    line.index(" <=")
+    for line in rendered.splitlines()
+    if " <=" in line
+] + [
+    line.index(" = ")
+    for line in rendered.splitlines()
+    if " = " in line
+]
+len(operator_columns) >= 2
+max(operator_columns) - min(operator_columns) <= 1
+```
+
 ```python doctest
->>> import arco
->>> model = arco.Model()
->>> x1 = model.add_variable(bounds=arco.Binary, name="x[1]")
->>> x2 = model.add_variable(bounds=arco.Binary, name="x[2]")
->>> x3 = model.add_variable(bounds=arco.Binary, name="x[3]")
->>> x4 = model.add_variable(bounds=arco.Binary, name="x[4]")
->>> x5 = model.add_variable(bounds=arco.Binary, name="x[5]")
->>> _ = model.add_constraint(2.0 * x1 + 8.0 * x2 + 4.0 * x3 + 2.0 * x4 + 5.0 * x5 <= 10.0)
->>> _ = model.add_constraint(x1 + x3 - x5 == 0.0)
->>> model.maximize(5.0 * x1 + 3.0 * x2 + 2.0 * x3 + 7.0 * x4 + 4.0 * x5)
->>> rendered = str(model)
->>> rendered.startswith("Max 5 x[1] + 3 x[2] + 2 x[3] + 7 x[4] + 4 x[5]")
-True
->>> "\ns.t.\n" in rendered
-True
->>> "Subject to" in rendered
-False
->>> "Binary: x[1], x[2], x[3], x[4], x[5]" in rendered
-True
->>> " <= 10" in rendered
-True
->>> operator_columns = [
-...     line.index(" <=")
-...     for line in rendered.splitlines()
-...     if " <=" in line
-... ] + [
-...     line.index(" = ")
-...     for line in rendered.splitlines()
-...     if " = " in line
-... ]
->>> len(operator_columns) >= 2
-True
->>> max(operator_columns) - min(operator_columns) <= 1
-True
+# >>> import arco
+# >>> model = arco.Model()
+# >>> x1 = model.add_variable(bounds=arco.Binary, name="x[1]")
+# >>> x2 = model.add_variable(bounds=arco.Binary, name="x[2]")
+# >>> x3 = model.add_variable(bounds=arco.Binary, name="x[3]")
+# >>> x4 = model.add_variable(bounds=arco.Binary, name="x[4]")
+# >>> x5 = model.add_variable(bounds=arco.Binary, name="x[5]")
+# >>> _ = model.add_constraint(2.0 * x1 + 8.0 * x2 + 4.0 * x3 + 2.0 * x4 + 5.0 * x5 <= 10.0)
+# >>> _ = model.add_constraint(x1 + x3 - x5 == 0.0)
+# >>> model.maximize(5.0 * x1 + 3.0 * x2 + 2.0 * x3 + 7.0 * x4 + 4.0 * x5)
+# >>> rendered = str(model)
+# >>> rendered.startswith("Max 5 x[1] + 3 x[2] + 2 x[3] + 7 x[4] + 4 x[5]")
+# True
+# >>> "\ns.t.\n" in rendered
+# True
+# >>> "Subject to" in rendered
+# False
+# >>> "Binary: x[1], x[2], x[3], x[4], x[5]" in rendered
+# True
+# >>> " <= 10" in rendered
+# True
+# >>> operator_columns = [
+# ...     line.index(" <=")
+# ...     for line in rendered.splitlines()
+# ...     if " <=" in line
+# ... ] + [
+# ...     line.index(" = ")
+# ...     for line in rendered.splitlines()
+# ...     if " = " in line
+# ... ]
+# >>> len(operator_columns) >= 2
+# True
+# >>> max(operator_columns) - min(operator_columns) <= 1
+# True
+```
+
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import contextlib
+import io
+import arco
+model = arco.Model()
+vars_ = [model.add_variable(bounds=arco.Binary, name=f"x[{i + 1}]") for i in range(35)]
+model.minimize(sum(vars_))
+for idx in range(22):
+    _ = model.add_constraint(sum(vars_) <= float(idx), name=f"c[{idx + 1}]")
+preview = str(model)
+"... (5 more terms)" in preview
+"... (2 more constraints)" in preview
+buf = io.StringIO()
+with contextlib.redirect_stdout(buf):
+    model.pprint()
+full = buf.getvalue()
+"... (2 more constraints)" in full
+"c[22]:" in full
 ```
 
 ```python doctest
->>> import contextlib
->>> import io
->>> import arco
->>> model = arco.Model()
->>> vars_ = [model.add_variable(bounds=arco.Binary, name=f"x[{i + 1}]") for i in range(35)]
->>> model.minimize(sum(vars_))
->>> for idx in range(22):
-...     _ = model.add_constraint(sum(vars_) <= float(idx), name=f"c[{idx + 1}]")
->>> preview = str(model)
->>> "... (5 more terms)" in preview
-True
->>> "... (2 more constraints)" in preview
-True
->>> buf = io.StringIO()
->>> with contextlib.redirect_stdout(buf):
-...     model.pprint()
->>> full = buf.getvalue()
->>> "... (2 more constraints)" in full
-False
->>> "c[22]:" in full
-True
+# >>> import contextlib
+# >>> import io
+# >>> import arco
+# >>> model = arco.Model()
+# >>> vars_ = [model.add_variable(bounds=arco.Binary, name=f"x[{i + 1}]") for i in range(35)]
+# >>> model.minimize(sum(vars_))
+# >>> for idx in range(22):
+# ...     _ = model.add_constraint(sum(vars_) <= float(idx), name=f"c[{idx + 1}]")
+# >>> preview = str(model)
+# >>> "... (5 more terms)" in preview
+# True
+# >>> "... (2 more constraints)" in preview
+# True
+# >>> buf = io.StringIO()
+# >>> with contextlib.redirect_stdout(buf):
+# ...     model.pprint()
+# >>> full = buf.getvalue()
+# >>> "... (2 more constraints)" in full
+# False
+# >>> "c[22]:" in full
+# True
+```
+
+<!-- copy-paste: auto-generated from doctest -->
+
+```python
+import arco
+model = arco.Model()
+t = arco.IndexSet(name="T", size=2)
+g = arco.IndexSet(name="G", members=["solar", "wind", "gas"])
+gen = model.add_variables(t, g, bounds=arco.Bounds(lower=0.0, upper=100.0), name="gen")
+capacity = {"solar": 50.0, "wind": 80.0, "gas": 100.0}
+caps = [capacity[name] for name in g.members] * t.size
+_ = model.add_constraints(gen <= caps)
+_ = model.add_constraints(gen.sum(over=g) >= [120.0, 90.0])
+rendered = str(model)
+"Index sets:" in rendered
+"T = [0, 1]" in rendered
+"G = [solar, wind, gas]" in rendered
+"gen[0,solar]" in rendered
+"gen[1,gas]" in rendered
+"Bounds:" in rendered
+"0 <= gen[t,g] <= 100  for t in T, g in G" in rendered
 ```
 
 ```python doctest
->>> import arco
->>> model = arco.Model()
->>> t = arco.IndexSet(name="T", size=2)
->>> g = arco.IndexSet(name="G", members=["solar", "wind", "gas"])
->>> gen = model.add_variables(t, g, bounds=arco.Bounds(lower=0.0, upper=100.0), name="gen")
->>> capacity = {"solar": 50.0, "wind": 80.0, "gas": 100.0}
->>> caps = [capacity[name] for name in g.members] * t.size
->>> _ = model.add_constraints(gen <= caps)
->>> _ = model.add_constraints(gen.sum(over=g) >= [120.0, 90.0])
->>> rendered = str(model)
->>> "Index sets:" in rendered
-True
->>> "T = [0, 1]" in rendered
-True
->>> "G = [solar, wind, gas]" in rendered
-True
->>> "gen[0,solar]" in rendered
-True
->>> "gen[1,gas]" in rendered
-True
->>> "Bounds:" in rendered
-True
->>> "0 <= gen[t,g] <= 100  for t in T, g in G" in rendered
-True
+# >>> import arco
+# >>> model = arco.Model()
+# >>> t = arco.IndexSet(name="T", size=2)
+# >>> g = arco.IndexSet(name="G", members=["solar", "wind", "gas"])
+# >>> gen = model.add_variables(t, g, bounds=arco.Bounds(lower=0.0, upper=100.0), name="gen")
+# >>> capacity = {"solar": 50.0, "wind": 80.0, "gas": 100.0}
+# >>> caps = [capacity[name] for name in g.members] * t.size
+# >>> _ = model.add_constraints(gen <= caps)
+# >>> _ = model.add_constraints(gen.sum(over=g) >= [120.0, 90.0])
+# >>> rendered = str(model)
+# >>> "Index sets:" in rendered
+# True
+# >>> "T = [0, 1]" in rendered
+# True
+# >>> "G = [solar, wind, gas]" in rendered
+# True
+# >>> "gen[0,solar]" in rendered
+# True
+# >>> "gen[1,gas]" in rendered
+# True
+# >>> "Bounds:" in rendered
+# True
+# >>> "0 <= gen[t,g] <= 100  for t in T, g in G" in rendered
+# True
 ```
diff --git a/scripts/test_docs_doctest.py b/scripts/test_docs_doctest.py
index e3ca063b..16b17e51 100644
--- a/scripts/test_docs_doctest.py
+++ b/scripts/test_docs_doctest.py
@@ -19,6 +19,8 @@ def _find_repo_root(*, start: Path) -> Path:
 REPO_ROOT = _find_repo_root(start=Path(__file__).resolve().parent)
 DOCS_DIR = REPO_ROOT / "docs"
 DOCTEST_FENCE_PATTERN = re.compile(r"^```python\s+doctest\b")
+COMMENTED_DOCTEST_PROMPT_PATTERN = re.compile(r"^\s*#\s*(>>>|\.\.\.)")
+COMMENTED_DOCTEST_LINE_PATTERN = re.compile(r"^(\s*)#\s?(.*)$")
 FENCE_END = "```"
 EXCLUDED_DOC_DIRS: set[str] = set()
 
@@ -94,14 +96,44 @@ def _collect_doctest_blocks() -> list[DoctestBlock]:
     return blocks
 
 
+def _normalize_doctest_source(*, source: str) -> str:
+    source_lines = source.splitlines()
+    uses_commented_doctest = any(
+        COMMENTED_DOCTEST_PROMPT_PATTERN.match(line) is not None
+        for line in source_lines
+        if line.strip()
+    )
+    if not uses_commented_doctest:
+        return source
+
+    normalized_lines: list[str] = []
+    for line_number, line in enumerate(source_lines, start=1):
+        if not line.strip():
+            normalized_lines.append(line)
+            continue
+
+        match = COMMENTED_DOCTEST_LINE_PATTERN.match(line)
+        if match is None:
+            raise AssertionError(
+                "Commented doctest block contains a non-comment line "
+                f"at source line {line_number}: {line!r}"
+            )
+
+        indent, content = match.groups()
+        normalized_lines.append(f"{indent}{content}")
+
+    return "\n".join(normalized_lines)
+
+
 DOCTEST_BLOCKS = _collect_doctest_blocks()
 
 
 @pytest.mark.parametrize("block", DOCTEST_BLOCKS, ids=lambda block: block.case_id)
 def test_markdown_doctest_blocks(block: DoctestBlock) -> None:
+    source = _normalize_doctest_source(source=block.source)
     parser = doctest.DocTestParser()
     document_test = parser.get_doctest(
-        block.source,
+        source,
         globs={"__name__": "__main__"},
         name=block.case_id,
         filename=str(block.file_path),
@@ -115,3 +147,26 @@ def test_markdown_doctest_blocks(block: DoctestBlock) -> None:
         f"Doctest failed for {block.case_id} "
         f"(attempted={result.attempted}, failed={result.failed})"
     )
+
+
+def test_normalize_doctest_source_keeps_standard_doctest() -> None:
+    source = ">>> answer = 2 + 2\n>>> answer\n4"
+
+    normalized = _normalize_doctest_source(source=source)
+
+    assert normalized == source
+
+
+def test_normalize_doctest_source_supports_comment_prefixed_blocks() -> None:
+    source = "# >>> answer = 2 + 2\n# >>> answer\n# 4"
+
+    normalized = _normalize_doctest_source(source=source)
+
+    assert normalized == ">>> answer = 2 + 2\n>>> answer\n4"
+
+
+def test_normalize_doctest_source_rejects_mixed_commented_blocks() -> None:
+    source = "# >>> answer = 2 + 2\n4"
+
+    with pytest.raises(AssertionError, match="non-comment line"):
+        _normalize_doctest_source(source=source)