"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Cooling duty of the second system is 62173345.96 kJ per year.\n",
"Heating duty of the second system is 54793314.40 kJ per year.\n"
]
}
],
"source": [
"# And many others related to heating and cooling utilities,\n",
"# since there is no use of heating/cooling duties in our system,\n",
"# let's just add two new units to demonstrate it\n",
"H1 = qs.unit_operations.HXutility('H1', ins=sys.outs[0], T=50+273.15)\n",
"H2 = qs.unit_operations.HXutility('H2', ins=H1-0, T=20+273.15)\n",
"sys2 = qs.System('sys2', path=(sys, H1, H2))\n",
"sys2.operating_hours = sys.operating_hours\n",
"sys2.simulate()\n",
"sys2.diagram()\n",
"print(f'Cooling duty of the second system is {sys2.get_cooling_duty():.2f} kJ per year.')\n",
"print(f'Heating duty of the second system is {sys2.get_heating_duty():.2f} kJ per year.')"
]
},
{
"cell_type": "markdown",
"id": "dc22ade4",
"metadata": {},
"source": [
"### 2.3. Inspecting and exporting results\n",
"Beyond the attributes above, a `System` offers a few high-value methods:\n",
"\n",
"- `sys.results()` returns a `DataFrame` summarizing every unit's design and cost.\n",
"- `sys.save_report('report.xlsx')` writes a full Excel report (stream tables, designs, costs, utilities) — the usual way to share results.\n",
"- `sys.diagram(kind='thorough', file='sys', format='png')` saves the flowsheet to a file; `kind` can be `'thorough'`, `'cluster'`, `'minimal'`, or `'surface'`."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "3f43d9fe",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:05.615878Z",
"iopub.status.busy": "2026-05-31T11:08:05.615878Z",
"iopub.status.idle": "2026-05-31T11:08:05.628710Z",
"shell.execute_reply": "2026-05-31T11:08:05.628710Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" System | \n",
" Units | \n",
" external_sys | \n",
"
\n",
" \n",
" \n",
" \n",
" | Electricity | \n",
" Power | \n",
" kW | \n",
" 0.836 | \n",
"
\n",
" \n",
" | Cost | \n",
" USD/hr | \n",
" 0.0654 | \n",
"
\n",
" \n",
" | Total purchase cost | \n",
" | \n",
" USD | \n",
" 8.61e+04 | \n",
"
\n",
" \n",
" | Installed equipment cost | \n",
" | \n",
" USD | \n",
" 1.42e+05 | \n",
"
\n",
" \n",
" | Utility cost | \n",
" | \n",
" USD/hr | \n",
" 0.0654 | \n",
"
\n",
" \n",
" | Material cost | \n",
" | \n",
" USD/hr | \n",
" 0 | \n",
"
\n",
" \n",
" | Sales | \n",
" | \n",
" USD/hr | \n",
" 0 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"System Units external_sys\n",
"Electricity Power kW 0.836\n",
" Cost USD/hr 0.0654\n",
"Total purchase cost USD 8.61e+04\n",
"Installed equipment cost USD 1.42e+05\n",
"Utility cost USD/hr 0.0654\n",
"Material cost USD/hr 0\n",
"Sales USD/hr 0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sys.results()"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "sys_export",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:05.631852Z",
"iopub.status.busy": "2026-05-31T11:08:05.630764Z",
"iopub.status.idle": "2026-05-31T11:08:05.635749Z",
"shell.execute_reply": "2026-05-31T11:08:05.634735Z"
}
},
"outputs": [],
"source": [
"# To export everything to a single Excel workbook (stream tables, unit designs,\n",
"# costs, and utilities), use `save_report`. It saves next to this notebook unless\n",
"# you pass a path via `file=`. The `results()` DataFrame can also be saved directly\n",
"# with pandas (e.g., `.to_csv(...)`). Commented out so no files are written here.\n",
"# sys.save_report('system_report.xlsx')\n",
"# sys.results().to_csv('system_results.csv')"
]
},
{
"cell_type": "markdown",
"id": "fbc02ff7",
"metadata": {},
"source": [
"\n",
"\n",
"**Heads up: dynamic simulation.** `QSDsan` also extends `System` to integrate unit states over time. With dynamic units you call `sys.simulate(t_span=(0, 10), ...)` instead of a steady-state solve. That is covered in [11. Dynamic Simulation](https://qsdsan.readthedocs.io/en/latest/tutorials/11_Dynamic_Simulation.html).\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "ed3f5fa1",
"metadata": {},
"source": [
"## 3. Process specifications \n",
"\n",
"Just like a unit (see the [SanUnit advanced tutorial section 3.6](5_SanUnit_advanced.ipynb#3.7.-Specifications)), a whole `System` can carry *specifications*: functions that run during `sys.simulate()`. Add one with `sys.add_specification`, often as a decorator, and pass `simulate=True` so the system re-converges after the spec has run. System-level specifications are the right altitude when the quantity you want to control depends on more than one unit, or on the converged recycle streams, rather than on a single unit in isolation."
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "2e7cdc77",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:05.637750Z",
"iopub.status.busy": "2026-05-31T11:08:05.637750Z",
"iopub.status.idle": "2026-05-31T11:08:05.717206Z",
"shell.execute_reply": "2026-05-31T11:08:05.716697Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Without a spec: product COD = 143.1 mg/L\n",
"With spec (target 200): product COD = 199.7 mg/L, dilution = 1150 kg/hr\n"
]
}
],
"source": [
"# Build a small blend-and-split train from the same kinds of units used above:\n",
"# a high-strength wastewater diluted with process water, then split to a product line.\n",
"conc = qs.WasteStream.codbased_inf_model('conc', flow_tot=1000, units=('L/hr', 'mg/L'))\n",
"dil = qs.WasteStream('dil', H2O=2000, units='kg/hr') # process water, negligible COD\n",
"M1 = qs.unit_operations.MixTank('M1', ins=(conc, dil), tau=1, V_wf=0.8, init_with='WasteStream')\n",
"S3 = qs.unit_operations.Splitter('S3', ins=M1-0, outs=('product', 'sidestream'),\n",
" split=0.7, init_with='WasteStream')\n",
"blend_sys = qs.System.from_units('blend_sys', units=[M1, S3])\n",
"blend_sys.simulate()\n",
"print(f'Without a spec: product COD = {qs.F.stream.product.COD:.1f} mg/L')\n",
"\n",
"# Hold the product COD at a target by adjusting the dilution-water flow.\n",
"# `simulate=True` re-converges the system after the spec sets the knob.\n",
"target_COD = 200.0 # mg/L\n",
"@blend_sys.add_specification(simulate=True)\n",
"def hold_product_COD():\n",
" needed_vol = conc.COD * conc.F_vol / target_COD # total volume to reach the target conc\n",
" dil.imass['H2O'] = max(needed_vol - conc.F_vol, 0.) * 1000 # ~1000 kg/m3 of added water\n",
"blend_sys.simulate()\n",
"print(f'With spec (target {target_COD:.0f}): product COD = {qs.F.stream.product.COD:.1f} mg/L, '\n",
" f'dilution = {dil.imass[\"H2O\"]:.0f} kg/hr')"
]
},
{
"cell_type": "markdown",
"id": "f3fd9e74",
"metadata": {},
"source": [
"Here one parameter (the dilution-water flow) is set to meet one target (the product COD). That balance is the whole story of a specification: **each spec consumes one degree of freedom.** To hit N independent targets you must free N parameters; asking for more targets than free parameter over-constrains the system and has no consistent solution. (This is the system-level version of the degrees-of-freedom point from the [SanUnit advanced tutorial section 3.6](5_SanUnit_advanced.ipynb#3.7.-Specifications).)\n",
"\n",
"We could compute the dilution analytically above because mixing is simple. In general an outlet depends on the converged system in ways you cannot invert by hand, and you instead drive a parameter to the target numerically. In this case, you can use `add_bounded_numerical_specification`, which adjusts one parameter until your objective function returns zero, using a bracketed root-finder."
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "b3626b1d",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:05.719220Z",
"iopub.status.busy": "2026-05-31T11:08:05.719220Z",
"iopub.status.idle": "2026-05-31T11:08:06.094128Z",
"shell.execute_reply": "2026-05-31T11:08:06.093110Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solved dilution = 1147 kg/hr -> product COD = 200.0 mg/L\n"
]
}
],
"source": [
"# Drive the same target by root-finding instead of a formula. First clear the analytic\n",
"# spec so the two do not both fight over the dilution knob (one knob, one spec).\n",
"blend_sys.specifications = []\n",
"\n",
"@blend_sys.add_bounded_numerical_specification(x0=0., x1=2e4, xtol=1.0)\n",
"def hit_product_COD(h2o):\n",
" dil.imass['H2O'] = h2o # the one knob: dilution-water flow (kg/hr)\n",
" blend_sys.simulate()\n",
" return qs.F.stream.product.COD - target_COD # objective: zero at the target\n",
"blend_sys.simulate()\n",
"print(f'Solved dilution = {dil.imass[\"H2O\"]:.0f} kg/hr -> product COD = {qs.F.stream.product.COD:.1f} mg/L')"
]
},
{
"cell_type": "markdown",
"id": "s4-intro",
"metadata": {},
"source": [
"## 4. Operational flexibility \n",
"\n",
"The systems built so far are each evaluated at a single steady state. Many real systems instead operate under conditions that vary over the year, such as seasonal loading, day and night cycles, or periodic feedstock switching. Evaluating a single representative condition can misstate the annual totals.\n",
"\n",
"An `AgileSystem` represents a system that runs in several operation modes over a year, and compiles one annualized result across all of them, with each mode weighted by its operating hours. It is re-exported from `biosteam` at the top level as `qs.AgileSystem`. This is particularly useful for techno-economic analysis (TEA) and life cycle assessment (LCA), where a single annualized result across operating modes is more representative than a single steady-state snapshot. Each operation mode is solved at steady state, so `AgileSystem` is not a dynamic (time-resolved) simulation; for transient behavior over time see [11. Dynamic Simulation](https://qsdsan.readthedocs.io/en/latest/tutorials/11_Dynamic_Simulation.html)."
]
},
{
"cell_type": "markdown",
"id": "s4-modes-md",
"metadata": {},
"source": [
"An agile system is assembled from two ingredients:\n",
"\n",
"- operation parameters: functions that set a quantity that varies between modes (for example, the influent strength), registered with `agile_sys.operation_parameter`; and\n",
"- operation modes: each a system paired with its operating hours and the parameter values that hold during that mode, registered with `agile_sys.operation_mode`.\n",
"\n",
"The individual mode objects are created through the `operation_mode` method rather than instantiated directly."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "s4-load",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:06.096146Z",
"iopub.status.busy": "2026-05-31T11:08:06.096146Z",
"iopub.status.idle": "2026-05-31T11:08:06.123190Z",
"shell.execute_reply": "2026-05-31T11:08:06.122182Z"
}
},
"outputs": [],
"source": [
"from qsdsan.utils import create_example_treatment_systems\n",
"\n",
"# Reuse the aerobic treatment plant from the TEA and LCA tutorials\n",
"# (municipal wastewater, 4,000 m3/d). A dedicated flowsheet isolates registries.\n",
"qs.main_flowsheet.set_flowsheet('agile_demo')\n",
"aer_sys, ana_sys = create_example_treatment_systems()\n",
"aer = aer_sys.units[0]\n",
"influent = aer.ins[0]\n",
"Substrate = influent.components.Substrate\n",
"design_flow = 4000/24 # m3/hr"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "s4-build",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:06.125197Z",
"iopub.status.busy": "2026-05-31T11:08:06.125197Z",
"iopub.status.idle": "2026-05-31T11:08:06.132476Z",
"shell.execute_reply": "2026-05-31T11:08:06.131456Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[OperationMode(system=aer_sys, operating_hours=4380.0, COD=800),\n",
" OperationMode(system=aer_sys, operating_hours=4380.0, COD=300)]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"agile_sys = qs.AgileSystem()\n",
"\n",
"@agile_sys.operation_parameter\n",
"def set_influent_COD(COD):\n",
" # influent strength in mg/L, applied as the substrate mass flow\n",
" influent.imass['Substrate'] = (COD/1000 * design_flow)/Substrate.i_COD\n",
"\n",
"half_year = 0.5 * 365 * 24 # hours\n",
"agile_sys.operation_mode(aer_sys, operating_hours=half_year, COD=800) # high-load season\n",
"agile_sys.operation_mode(aer_sys, operating_hours=half_year, COD=300) # low-load season\n",
"\n",
"agile_sys.simulate()\n",
"agile_sys.operation_modes"
]
},
{
"cell_type": "markdown",
"id": "s4-results-md",
"metadata": {},
"source": [
"After simulation, the agile system reports results aggregated across all modes. The aerobic plant spends electricity on aeration in proportion to the organic load it oxidizes, so the higher-load season draws more power. The annual electricity consumption is the sum over modes of each mode's power multiplied by its operating hours, which we confirm below."
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "s4-results",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:06.134481Z",
"iopub.status.busy": "2026-05-31T11:08:06.134481Z",
"iopub.status.idle": "2026-05-31T11:08:06.140474Z",
"shell.execute_reply": "2026-05-31T11:08:06.139351Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Total operating hours: 8760 hr\n",
"Annual electricity consumption: 381425 kWh\n",
" COD 800 mg/L: 63.3 kW over 4380 hr\n",
" COD 300 mg/L: 23.7 kW over 4380 hr\n"
]
}
],
"source": [
"print(f'Total operating hours: {agile_sys.operating_hours:.0f} hr')\n",
"print(f'Annual electricity consumption: {agile_sys.get_electricity_consumption():.0f} kWh')\n",
"\n",
"# The annualized value is the operating-hour-weighted sum over the modes\n",
"for mode in agile_sys.operation_modes:\n",
" set_influent_COD(mode.COD); aer_sys.simulate()\n",
" print(f' COD {mode.COD:>4} mg/L: {aer.power_utility.rate:5.1f} kW '\n",
" f'over {mode.operating_hours:.0f} hr')\n",
"agile_sys.simulate() # restore the aggregated state"
]
},
{
"cell_type": "markdown",
"id": "s4-tea-md",
"metadata": {},
"source": [
"The TEA and LCA classes accept an agile system in place of a single system, so one analysis spans all modes. Operating costs that depend on the mode, such as electricity, reflect the operating-hour-weighted result, while the installed capital is shared across modes."
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "s4-tea",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:06.143482Z",
"iopub.status.busy": "2026-05-31T11:08:06.142481Z",
"iopub.status.idle": "2026-05-31T11:08:06.331444Z",
"shell.execute_reply": "2026-05-31T11:08:06.330422Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Installed equipment cost: 5,000,000 USD\n",
"Annual operating cost: 48,145 USD/yr\n",
"Net present value: -10,599,997 USD\n"
]
}
],
"source": [
"tea = qs.TEA(agile_sys, discount_rate=0.05, lifetime=20, simulate_system=False)\n",
"print(f'Installed equipment cost: {agile_sys.installed_equipment_cost:,.0f} USD')\n",
"print(f'Annual operating cost: {tea.AOC:,.0f} USD/yr')\n",
"print(f'Net present value: {tea.NPV:,.0f} USD')"
]
},
{
"cell_type": "markdown",
"id": "s4-lca-md",
"metadata": {},
"source": [
"An LCA spans all modes in the same way. Construction impacts are one-time and independent of the operating mode, whereas operational impacts such as electricity use the annualized consumption reported by the agile system. Below we attach construction materials to the plant and supply the annualized electricity, each characterized for global warming potential."
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "s4-lca",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:06.333443Z",
"iopub.status.busy": "2026-05-31T11:08:06.333443Z",
"iopub.status.idle": "2026-05-31T11:08:06.341108Z",
"shell.execute_reply": "2026-05-31T11:08:06.339910Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Total GWP over 20 yr: 3,964,250 kg CO2-eq\n",
" construction: 150,000\n",
" operation: 3,814,250\n",
"Annualized GWP: 198,212 kg CO2-eq/yr\n"
]
}
],
"source": [
"GWP = qs.ImpactIndicator('GlobalWarming', alias='GWP', unit='kg CO2-eq')\n",
"Concrete = qs.ImpactItem('Concrete', functional_unit='m3', GWP=300.)\n",
"Steel = qs.ImpactItem('Steel', functional_unit='kg', GWP=2.)\n",
"electricity = qs.ImpactItem('electricity', functional_unit='kWh', GWP=0.5)\n",
"\n",
"aer.construction = [\n",
" qs.Construction(linked_unit=aer, item=Concrete, quantity=300, quantity_unit='m3', lifetime=30),\n",
" qs.Construction(linked_unit=aer, item=Steel, quantity=30000, quantity_unit='kg', lifetime=30),\n",
"]\n",
"\n",
"lifetime = 20\n",
"annual_electricity = agile_sys.get_electricity_consumption() # kWh/yr, across all modes\n",
"lca = qs.LCA(system=agile_sys, lifetime=lifetime,\n",
" electricity=annual_electricity*lifetime, simulate_system=False)\n",
"print(f'Total GWP over {lifetime} yr: {lca.get_total_impacts()[\"GlobalWarming\"]:,.0f} kg CO2-eq')\n",
"print(f' construction: {lca.total_construction_impacts[\"GlobalWarming\"]:,.0f}')\n",
"print(f' operation: {lca.get_total_impacts(operation_only=True)[\"GlobalWarming\"]:,.0f}')\n",
"print(f'Annualized GWP: {lca.get_total_impacts(annual=True)[\"GlobalWarming\"]:,.0f} kg CO2-eq/yr')"
]
},
{
"cell_type": "markdown",
"id": "s4-pointers",
"metadata": {},
"source": [
"\n",
"\n",
"**Mode-dependent parameters and further analysis.** An operation parameter can also depend on the mode itself (pass `mode_dependent=True` to `operation_parameter`), and metrics can be summed across modes with `annual_operation_metrics`. See `?qs.AgileSystem` for the full interface. For how the annualized TEA and LCA results are read, interpreted, and exported, see [Tutorial 7 on techno-economic analysis](7_TEA.ipynb) and [Tutorial 8 on life cycle assessment](8_LCA.ipynb).\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "conv-s5",
"metadata": {},
"source": [
"## 5. Controlling recycle convergence \n",
"\n",
"A system with recycle streams is solved iteratively: each unit runs in order, the recycle streams are updated, and the loop repeats until the recycle streams stop changing (within a tolerance). The systems above converged on their own, which is the common case. This section gives some tips on what to do when one does not: how to read the convergence report, the settings that control the solve and their defaults, and what to change first. It is the convergence-control surface that the recycle-convergence note in [14. Modeling Notes & Pitfalls](14_Modeling_Notes_and_Pitfalls.ipynb) points to.\n",
"\n",
"\n",
"\n",
"**Two different kinds of convergence.** This section is about steady-state recycle convergence, the iteration that resolves recycle loops in `simulate()`. Dynamic systems integrate unit states over time instead, where convergence means the ODE solver advancing through the time span. Those settings (the integration method, the time-step tolerances) are separate and covered in [11. Dynamic Simulation](https://qsdsan.readthedocs.io/en/latest/tutorials/11_Dynamic_Simulation.html).\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "conv-5_1-md",
"metadata": {},
"source": [
"### 5.1. Reading the convergence report\n",
"\n",
"After each `simulate()`, the system compares the recycle streams against the tolerances and reports the largest remaining error. This is the block shown under the diagram in Section 2: the first line is the worst component flow-rate error (in kmol/hr, with its value as a percentage of the stream flow in parentheses), the second is the temperature error (in K, with its percentage), and \"after N loops\" is the number of iterations taken.\n",
"\n",
"When the solver reaches `maxiter` without meeting the tolerances, `simulate()` raises a `RuntimeError` carrying that same report, prefixed with \"could not converge\". The cell below forces that error on purpose with far too strict tolerances and too few iterations, so the traceback below is the error. The cell after it restores the defaults and converges normally again."
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "conv-5_1-code",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:06.343118Z",
"iopub.status.busy": "2026-05-31T11:08:06.343118Z",
"iopub.status.idle": "2026-05-31T11:08:07.516727Z",
"shell.execute_reply": "2026-05-31T11:08:07.516727Z"
},
"tags": [
"raises-exception"
]
},
"outputs": [
{
"ename": "RuntimeError",
"evalue": " could not converge\nHighest convergence error among components in recycle\nstream S2-1 after 2 loops:\n- flow rate 2.04e-01 kmol/hr (20%)\n- temperature 0.00e+00 K (0%)",
"output_type": "error",
"traceback": [
"\u001b[31m---------------------------------------------------------------------------\u001b[39m",
"\u001b[31mRuntimeError\u001b[39m Traceback (most recent call last)",
"\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[30]\u001b[39m\u001b[32m, line 10\u001b[39m\n\u001b[32m 6\u001b[39m \n\u001b[32m 7\u001b[39m sys.molar_tolerance = sys.relative_molar_tolerance = \u001b[32m1e-9\u001b[39m \u001b[38;5;66;03m# far too strict\u001b[39;00m\n\u001b[32m 8\u001b[39m sys.temperature_tolerance = sys.relative_temperature_tolerance = \u001b[32m1e-9\u001b[39m\n\u001b[32m 9\u001b[39m sys.maxiter = \u001b[32m2\u001b[39m \u001b[38;5;66;03m# far too few iterations\u001b[39;00m\n\u001b[32m---> \u001b[39m\u001b[32m10\u001b[39m sys.simulate()\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\biosteam\\_system.py:3374\u001b[39m, in \u001b[36mSystem.simulate\u001b[39m\u001b[34m(self, update_configuration, units, design_and_cost, **kwargs)\u001b[39m\n\u001b[32m 3354\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34msimulate\u001b[39m(\u001b[38;5;28mself\u001b[39m, update_configuration: Optional[\u001b[38;5;28mbool\u001b[39m]=\u001b[38;5;28;01mNone\u001b[39;00m, units=\u001b[38;5;28;01mNone\u001b[39;00m, \n\u001b[32m 3355\u001b[39m design_and_cost=\u001b[38;5;28;01mNone\u001b[39;00m, **kwargs):\n\u001b[32m 3356\u001b[39m \u001b[38;5;250m \u001b[39m\u001b[33;03m\"\"\"\u001b[39;00m\n\u001b[32m 3357\u001b[39m \u001b[33;03m If system is dynamic, run the system dynamically. Otherwise, converge \u001b[39;00m\n\u001b[32m 3358\u001b[39m \u001b[33;03m the path of unit operations to steady state. After running/converging \u001b[39;00m\n\u001b[32m (...)\u001b[39m\u001b[32m 3372\u001b[39m \u001b[33;03m \u001b[39;00m\n\u001b[32m 3373\u001b[39m \u001b[33;03m \"\"\"\u001b[39;00m\n\u001b[32m-> \u001b[39m\u001b[32m3374\u001b[39m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m.flowsheet:\n\u001b[32m 3375\u001b[39m specifications = \u001b[38;5;28mself\u001b[39m._specifications\n\u001b[32m 3376\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m specifications \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m._running_specifications:\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\biosteam\\_flowsheet.py:120\u001b[39m, in \u001b[36mFlowsheet.__exit__\u001b[39m\u001b[34m(self, type, exception, traceback)\u001b[39m\n\u001b[32m 118\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m__exit__\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;28mtype\u001b[39m, exception, traceback):\n\u001b[32m 119\u001b[39m main_flowsheet.set_flowsheet(\u001b[38;5;28mself\u001b[39m._temporary_stack.pop())\n\u001b[32m--> \u001b[39m\u001b[32m120\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m exception: \u001b[38;5;28;01mraise\u001b[39;00m exception\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\biosteam\\_system.py:3435\u001b[39m, in \u001b[36mSystem.simulate\u001b[39m\u001b[34m(self, update_configuration, units, design_and_cost, **kwargs)\u001b[39m\n\u001b[32m 3429\u001b[39m outputs = \u001b[38;5;28mself\u001b[39m.simulate(\n\u001b[32m 3430\u001b[39m update_configuration=\u001b[38;5;28;01mTrue\u001b[39;00m, \n\u001b[32m 3431\u001b[39m design_and_cost=design_and_cost,\n\u001b[32m 3432\u001b[39m **kwargs\n\u001b[32m 3433\u001b[39m )\n\u001b[32m 3434\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m-> \u001b[39m\u001b[32m3435\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m error\n\u001b[32m 3436\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m 3437\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m (\u001b[38;5;129;01mnot\u001b[39;00m update_configuration \u001b[38;5;66;03m# Avoid infinite loop\u001b[39;00m\n\u001b[32m 3438\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mself\u001b[39m._connections != \u001b[38;5;28mself\u001b[39m._get_connections()):\n\u001b[32m 3439\u001b[39m \u001b[38;5;66;03m# Connections has been updated within simulation.\u001b[39;00m\n\u001b[32m 3440\u001b[39m \u001b[38;5;66;03m# Must reset path.\u001b[39;00m\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\biosteam\\_system.py:3422\u001b[39m, in \u001b[36mSystem.simulate\u001b[39m\u001b[34m(self, update_configuration, units, design_and_cost, **kwargs)\u001b[39m\n\u001b[32m 3420\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m 3421\u001b[39m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[32m-> \u001b[39m\u001b[32m3422\u001b[39m outputs = \u001b[30;43mself\u001b[39;49m\u001b[30;43m.\u001b[39;49m\u001b[30;43mconverge\u001b[39;49m\u001b[30;43m(\u001b[39;49m\u001b[30;43m*\u001b[39;49m\u001b[30;43m*\u001b[39;49m\u001b[30;43mkwargs\u001b[39;49m\u001b[30;43m)\u001b[39;49m\n\u001b[32m 3423\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m design_and_cost: \u001b[38;5;28mself\u001b[39m._summary()\n\u001b[32m 3424\u001b[39m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m error:\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\biosteam\\_system.py:3046\u001b[39m, in \u001b[36mSystem.converge\u001b[39m\u001b[34m(self, recycle_data, update_recycle_data)\u001b[39m\n\u001b[32m 3044\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;28mself\u001b[39m._N_runs): method()\n\u001b[32m 3045\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m-> \u001b[39m\u001b[32m3046\u001b[39m \u001b[30;43mmethod\u001b[39;49m\u001b[30;43m(\u001b[39;49m\u001b[30;43m)\u001b[39;49m\n\u001b[32m 3047\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m update_recycle_data:\n\u001b[32m 3048\u001b[39m \u001b[38;5;28;01mtry\u001b[39;00m: recycle_data.update()\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\biosteam\\_system.py:3002\u001b[39m, in \u001b[36mSystem._solve\u001b[39m\u001b[34m(self)\u001b[39m\n\u001b[32m 3000\u001b[39m data = \u001b[38;5;28mself\u001b[39m._get_recycle_data()\n\u001b[32m 3001\u001b[39m f = \u001b[38;5;28mself\u001b[39m._iter_run_conditional \u001b[38;5;28;01mif\u001b[39;00m conditional \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;28mself\u001b[39m._iter_run\n\u001b[32m-> \u001b[39m\u001b[32m3002\u001b[39m \u001b[38;5;28;01mtry\u001b[39;00m: \u001b[30;43msolver\u001b[39;49m\u001b[30;43m(\u001b[39;49m\u001b[30;43mf\u001b[39;49m\u001b[30;43m,\u001b[39;49m\u001b[30;43m \u001b[39;49m\u001b[30;43mdata\u001b[39;49m\u001b[30;43m,\u001b[39;49m\u001b[30;43m \u001b[39;49m\u001b[30;43m*\u001b[39;49m\u001b[30;43m*\u001b[39;49m\u001b[30;43mkwargs\u001b[39;49m\u001b[30;43m)\u001b[39;49m\n\u001b[32m 3003\u001b[39m \u001b[38;5;28;01mexcept\u001b[39;00m (\u001b[38;5;167;01mIndexError\u001b[39;00m, \u001b[38;5;167;01mValueError\u001b[39;00m) \u001b[38;5;28;01mas\u001b[39;00m error:\n\u001b[32m 3004\u001b[39m data = \u001b[38;5;28mself\u001b[39m._get_recycle_data()\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\flexsolve\\fixed_point_solvers.py:188\u001b[39m, in \u001b[36mconditional_aitken\u001b[39m\u001b[34m(f, x)\u001b[39m\n\u001b[32m 186\u001b[39m g, condition = f(x)\n\u001b[32m 187\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m condition: \u001b[38;5;28;01mreturn\u001b[39;00m g\n\u001b[32m--> \u001b[39m\u001b[32m188\u001b[39m gg, condition = \u001b[30;43mf\u001b[39;49m\u001b[30;43m(\u001b[39;49m\u001b[30;43mg\u001b[39;49m\u001b[30;43m)\u001b[39;49m\n\u001b[32m 189\u001b[39m x = aitken_iter(x, gg, x - g, gg - g)\n\u001b[32m 190\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m x\n",
"\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\.venv\\Lib\\site-packages\\biosteam\\_system.py:2403\u001b[39m, in \u001b[36mSystem._iter_run_conditional\u001b[39m\u001b[34m(self, data)\u001b[39m\n\u001b[32m 2400\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m (\u001b[38;5;28mself\u001b[39m._iter >= \u001b[38;5;28mself\u001b[39m.maxiter\n\u001b[32m 2401\u001b[39m \u001b[38;5;129;01mor\u001b[39;00m (\u001b[38;5;28mself\u001b[39m.maxtime \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mself\u001b[39m.timer.elapsed_time > \u001b[38;5;28mself\u001b[39m.maxtime)): \n\u001b[32m 2402\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m.strict_convergence: \n\u001b[32m-> \u001b[39m\u001b[32m2403\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mRuntimeError\u001b[39;00m(\u001b[33mf\u001b[39m\u001b[33m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mrepr\u001b[39m(\u001b[38;5;28mself\u001b[39m)\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m could not converge\u001b[39m\u001b[33m'\u001b[39m + \u001b[38;5;28mself\u001b[39m._error_info())\n\u001b[32m 2404\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m: \n\u001b[32m 2405\u001b[39m not_converged = \u001b[38;5;28;01mFalse\u001b[39;00m\n",
"\u001b[31mRuntimeError\u001b[39m: could not converge\nHighest convergence error among components in recycle\nstream S2-1 after 2 loops:\n- flow rate 2.04e-01 kmol/hr (20%)\n- temperature 0.00e+00 K (0%)"
]
}
],
"source": [
"# Reuse the system built above (sys = external_sys). Save the defaults so the next\n",
"# cell can restore them, then force a non-convergence error on purpose: with far too\n",
"# strict tolerances and too few iterations, simulate() cannot converge and raises.\n",
"saved = (sys.molar_tolerance, sys.relative_molar_tolerance,\n",
" sys.temperature_tolerance, sys.relative_temperature_tolerance, sys.maxiter)\n",
"\n",
"sys.molar_tolerance = sys.relative_molar_tolerance = 1e-9 # far too strict\n",
"sys.temperature_tolerance = sys.relative_temperature_tolerance = 1e-9\n",
"sys.maxiter = 2 # far too few iterations\n",
"sys.simulate()"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "conv-5_1-restore",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:07.516727Z",
"iopub.status.busy": "2026-05-31T11:08:07.516727Z",
"iopub.status.idle": "2026-05-31T11:08:07.525928Z",
"shell.execute_reply": "2026-05-31T11:08:07.525928Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Re-converged with the default settings.\n"
]
}
],
"source": [
"# Restore the saved settings and converge normally again.\n",
"(sys.molar_tolerance, sys.relative_molar_tolerance,\n",
" sys.temperature_tolerance, sys.relative_temperature_tolerance, sys.maxiter) = saved\n",
"sys.simulate()\n",
"print('Re-converged with the default settings.')"
]
},
{
"cell_type": "markdown",
"id": "conv-5_2-md",
"metadata": {},
"source": [
"### 5.2. The settings that control the solve\n",
"\n",
"Recycle convergence is governed by a handful of `System` attributes:\n",
"\n",
"- `molar_tolerance` (kmol/hr) and `relative_molar_tolerance` (a fraction): the absolute and relative flow-rate tolerances.\n",
"- `temperature_tolerance` (K) and `relative_temperature_tolerance` (a fraction): the absolute and relative temperature tolerances.\n",
"- `maxiter`: the maximum number of iterations before `simulate()` gives up and raises exceptions.\n",
"- `method` (also exposed as `converge_method`): the acceleration method used to update the recycle guess between iterations.\n",
"\n",
"The solve stops once a recycle is within the molar tolerance (absolute or relative) and within the temperature tolerance (absolute or relative). Set any of these attributes directly, or set the common ones together with `set_tolerance(mol=, rmol=, T=, maxiter=)`. The live defaults and the available methods are printed below."
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "conv-5_2-code",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:07.528256Z",
"iopub.status.busy": "2026-05-31T11:08:07.528256Z",
"iopub.status.idle": "2026-05-31T11:08:07.534828Z",
"shell.execute_reply": "2026-05-31T11:08:07.534828Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"molar_tolerance = 1.0 kmol/hr\n",
"relative_molar_tolerance = 0.01\n",
"temperature_tolerance = 0.1 K\n",
"relative_temperature_tolerance = 0.001\n",
"maxiter = 200\n",
"method = 'aitken'\n",
"available methods = ['aitken', 'wegstein', 'fixedpoint', 'anderson', 'diagbroyden', 'excitingmixing', 'linearmixing', 'broyden1', 'broyden2']\n"
]
}
],
"source": [
"print(f'molar_tolerance = {sys.molar_tolerance} kmol/hr')\n",
"print(f'relative_molar_tolerance = {sys.relative_molar_tolerance}')\n",
"print(f'temperature_tolerance = {sys.temperature_tolerance} K')\n",
"print(f'relative_temperature_tolerance = {sys.relative_temperature_tolerance}')\n",
"print(f'maxiter = {sys.maxiter}')\n",
"print(f'method = {sys.method!r}')\n",
"print(f'available methods = {list(sys.available_methods)}')"
]
},
{
"cell_type": "markdown",
"id": "conv-5_3-md",
"metadata": {},
"source": [
"### 5.3. What to adjust when a system will not converge\n",
"\n",
"When a system stalls or raises a convergence error, change things roughly in this order:\n",
"\n",
"1. **Loosen the tolerances** if the defaults are tighter than the analysis needs. The default `molar_tolerance` of 1.0 kmol/hr and `temperature_tolerance` of 0.1 K may be stricter than many wastewater analyses require; a larger value is looser and lets the solve stop sooner.\n",
"2. **Improve the recycle initial guess.** The solver starts from whatever is already on the recycle streams, so writing realistic flows onto `sys.recycle` before `simulate()` gives it a closer starting point. This is often the single most effective change for a stiff biokinetic loop.\n",
"3. **Raise `maxiter`** if the loop is converging but has not finished in time.\n",
"4. **Switch the acceleration method.** The default `'aitken'` is fastest near the solution but can over-correct when the guess is far off; `'wegstein'` is often steadier for an oscillating loop, and `'fixedpoint'` is plain successive substitution (no acceleration), the most robust but slowest."
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "conv-5_3-code",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-31T11:08:07.534828Z",
"iopub.status.busy": "2026-05-31T11:08:07.534828Z",
"iopub.status.idle": "2026-05-31T11:08:07.611919Z",
"shell.execute_reply": "2026-05-31T11:08:07.611919Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Highest convergence error among components in recycle\n",
"stream S2-1 after 2 loops:\n",
"- flow rate 4.82e+00 kmol/hr (39%)\n",
"- temperature 0.00e+00 K (0%)\n"
]
}
],
"source": [
"# 1. Loosen the tolerances (a larger number is looser than the default).\n",
"sys.set_tolerance(mol=5.0, T=0.5) # defaults are 1.0 kmol/hr and 0.1 K\n",
"\n",
"# 2. Seed the recycle stream with a realistic guess before solving.\n",
"sys.recycle.imass['H2O'] = ww.imass['H2O'] # the clarifier recycle (S2-1)\n",
"\n",
"# 3. Allow more iterations if needed.\n",
"sys.maxiter = 400\n",
"\n",
"# 4. Try a steadier acceleration method.\n",
"sys.method = 'wegstein'\n",
"\n",
"sys.simulate()\n",
"print(sys._error_info())"
]
},
{
"cell_type": "markdown",
"id": "nav-footer-6_system",
"metadata": {},
"source": [
"\n",
"\n",
"---\n",
"\n",
"↑ Back to top\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.7"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 5
}