{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using the ``ws3`` optimization modelling function (advanced) \n", "\n", "The `ws3.forest.ForestMdoel` class includes a number of methods to automate formulation and solution of optimization problems. In this notebook, we will demonstrate how to use these methods to formulate various common forms of objective functions and constraints, how to combine these in various ways, and how to solve these problems using the built-in solvers." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set up environment" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below we (optionally) install the `ws3` package from the local source code in this repository using the `-e` flag to ensure that the package is installed in editable mode (i.e., any changes you make to the source code immediately affect `ws3` behaviour the next time you run the notebook). This is is not necessary if you have installed `ws3` using pip or another method.\n", "\n", "Note that the code in this notebook has been tested with the version of the `ws3` package that is included in this repository. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set auto-reload to reload modules when they are changed." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Found existing installation: ws3 1.1.0.dev0\n", "Uninstalling ws3-1.1.0.dev0:\n", " Successfully uninstalled ws3-1.1.0.dev0\n", "Note: you may need to restart the kernel to use updated packages.\n", "Obtaining file:///home/gep/tmp/ws3\n", " Installing build dependencies ... \u001b[?25ldone\n", "\u001b[?25h Checking if build backend supports build_editable ... \u001b[?25ldone\n", "\u001b[?25h Getting requirements to build editable ... \u001b[?25ldone\n", "\u001b[?25h Installing backend dependencies ... \u001b[?25ldone\n", "\u001b[?25h Preparing editable metadata (pyproject.toml) ... \u001b[?25ldone\n", "\u001b[?25hRequirement already satisfied: dill in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (0.4.0)\n", "Requirement already satisfied: fiona in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (1.10.1)\n", "Requirement already satisfied: gurobipy in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (12.0.3)\n", "Requirement already satisfied: highspy in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (1.11.0)\n", "Requirement already satisfied: libcbm in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (2.8.1)\n", "Requirement already satisfied: matplotlib in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (3.10.5)\n", "Requirement already satisfied: numpy>=1.21 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (2.2.6)\n", "Requirement already satisfied: pandas>=1.3 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (2.3.1)\n", "Requirement already satisfied: profilehooks in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (1.13.0)\n", "Requirement already satisfied: pulp in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (3.2.2)\n", "Requirement already satisfied: rasterio in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (1.4.3)\n", "Requirement already satisfied: scipy>=1.7 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ws3==1.1.0.dev0) (1.16.1)\n", "Requirement already satisfied: python-dateutil>=2.8.2 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from pandas>=1.3->ws3==1.1.0.dev0) (2.9.0.post0)\n", "Requirement already satisfied: pytz>=2020.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from pandas>=1.3->ws3==1.1.0.dev0) (2025.2)\n", "Requirement already satisfied: tzdata>=2022.7 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from pandas>=1.3->ws3==1.1.0.dev0) (2025.2)\n", "Requirement already satisfied: six>=1.5 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from python-dateutil>=2.8.2->pandas>=1.3->ws3==1.1.0.dev0) (1.17.0)\n", "Requirement already satisfied: attrs>=19.2.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from fiona->ws3==1.1.0.dev0) (25.3.0)\n", "Requirement already satisfied: certifi in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from fiona->ws3==1.1.0.dev0) (2025.8.3)\n", "Requirement already satisfied: click~=8.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from fiona->ws3==1.1.0.dev0) (8.2.1)\n", "Requirement already satisfied: click-plugins>=1.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from fiona->ws3==1.1.0.dev0) (1.1.1.2)\n", "Requirement already satisfied: cligj>=0.5 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from fiona->ws3==1.1.0.dev0) (0.7.2)\n", "Requirement already satisfied: numexpr>=2.8.7 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from libcbm->ws3==1.1.0.dev0) (2.11.0)\n", "Requirement already satisfied: numba in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from libcbm->ws3==1.1.0.dev0) (0.61.2)\n", "Requirement already satisfied: pyyaml in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from libcbm->ws3==1.1.0.dev0) (6.0.2)\n", "Requirement already satisfied: mock in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from libcbm->ws3==1.1.0.dev0) (5.2.0)\n", "Requirement already satisfied: openpyxl in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from libcbm->ws3==1.1.0.dev0) (3.1.5)\n", "Requirement already satisfied: contourpy>=1.0.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib->ws3==1.1.0.dev0) (1.3.3)\n", "Requirement already satisfied: cycler>=0.10 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib->ws3==1.1.0.dev0) (0.12.1)\n", "Requirement already satisfied: fonttools>=4.22.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib->ws3==1.1.0.dev0) (4.59.0)\n", "Requirement already satisfied: kiwisolver>=1.3.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib->ws3==1.1.0.dev0) (1.4.8)\n", "Requirement already satisfied: packaging>=20.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib->ws3==1.1.0.dev0) (25.0)\n", "Requirement already satisfied: pillow>=8 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib->ws3==1.1.0.dev0) (11.3.0)\n", "Requirement already satisfied: pyparsing>=2.3.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib->ws3==1.1.0.dev0) (3.2.3)\n", "Requirement already satisfied: llvmlite<0.45,>=0.44.0dev0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from numba->libcbm->ws3==1.1.0.dev0) (0.44.0)\n", "Requirement already satisfied: et-xmlfile in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from openpyxl->libcbm->ws3==1.1.0.dev0) (2.0.0)\n", "Requirement already satisfied: affine in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from rasterio->ws3==1.1.0.dev0) (2.4.0)\n", "Building wheels for collected packages: ws3\n", " Building editable for ws3 (pyproject.toml) ... \u001b[?25ldone\n", "\u001b[?25h Created wheel for ws3: filename=ws3-1.1.0.dev0-py3-none-any.whl size=4136 sha256=5df49d1809d2f15004f3046010971c142390fb47fb3bfa7620b09e11903424a7\n", " Stored in directory: /tmp/pip-ephem-wheel-cache-fjqbp0j3/wheels/fa/4b/10/3fe4b92a02fb87987a6fe53a10fad0a22a781bf98cd7b63f17\n", "Successfully built ws3\n", "Installing collected packages: ws3\n", "Successfully installed ws3-1.1.0.dev0\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "clobber_ws3 = True\n", "if clobber_ws3:\n", " %pip uninstall -y ws3\n", " %pip install -e .." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use `pip` to install Python packages listed in `requirements.txt` (some extra packages needed for example notebooks to run correctly)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: seaborn in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from -r requirements.txt (line 1)) (0.13.2)\n", "Requirement already satisfied: geopandas in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from -r requirements.txt (line 2)) (1.1.1)\n", "Requirement already satisfied: ipywidgets in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from -r requirements.txt (line 3)) (8.1.7)\n", "Requirement already satisfied: numpy!=1.24.0,>=1.20 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from seaborn->-r requirements.txt (line 1)) (2.2.6)\n", "Requirement already satisfied: pandas>=1.2 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from seaborn->-r requirements.txt (line 1)) (2.3.1)\n", "Requirement already satisfied: matplotlib!=3.6.1,>=3.4 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from seaborn->-r requirements.txt (line 1)) (3.10.5)\n", "Requirement already satisfied: pyogrio>=0.7.2 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from geopandas->-r requirements.txt (line 2)) (0.11.1)\n", "Requirement already satisfied: packaging in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from geopandas->-r requirements.txt (line 2)) (25.0)\n", "Requirement already satisfied: pyproj>=3.5.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from geopandas->-r requirements.txt (line 2)) (3.7.1)\n", "Requirement already satisfied: shapely>=2.0.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from geopandas->-r requirements.txt (line 2)) (2.1.1)\n", "Requirement already satisfied: comm>=0.1.3 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipywidgets->-r requirements.txt (line 3)) (0.2.3)\n", "Requirement already satisfied: ipython>=6.1.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipywidgets->-r requirements.txt (line 3)) (9.4.0)\n", "Requirement already satisfied: traitlets>=4.3.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipywidgets->-r requirements.txt (line 3)) (5.14.3)\n", "Requirement already satisfied: widgetsnbextension~=4.0.14 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipywidgets->-r requirements.txt (line 3)) (4.0.14)\n", "Requirement already satisfied: jupyterlab_widgets~=3.0.15 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipywidgets->-r requirements.txt (line 3)) (3.0.15)\n", "Requirement already satisfied: decorator in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (5.2.1)\n", "Requirement already satisfied: ipython-pygments-lexers in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (1.1.1)\n", "Requirement already satisfied: jedi>=0.16 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (0.19.2)\n", "Requirement already satisfied: matplotlib-inline in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (0.1.7)\n", "Requirement already satisfied: pexpect>4.3 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (4.9.0)\n", "Requirement already satisfied: prompt_toolkit<3.1.0,>=3.0.41 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (3.0.51)\n", "Requirement already satisfied: pygments>=2.4.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (2.19.2)\n", "Requirement already satisfied: stack_data in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (0.6.3)\n", "Requirement already satisfied: wcwidth in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from prompt_toolkit<3.1.0,>=3.0.41->ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (0.2.13)\n", "Requirement already satisfied: parso<0.9.0,>=0.8.4 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from jedi>=0.16->ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (0.8.4)\n", "Requirement already satisfied: contourpy>=1.0.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (1.3.3)\n", "Requirement already satisfied: cycler>=0.10 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (0.12.1)\n", "Requirement already satisfied: fonttools>=4.22.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (4.59.0)\n", "Requirement already satisfied: kiwisolver>=1.3.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (1.4.8)\n", "Requirement already satisfied: pillow>=8 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (11.3.0)\n", "Requirement already satisfied: pyparsing>=2.3.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (3.2.3)\n", "Requirement already satisfied: python-dateutil>=2.7 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (2.9.0.post0)\n", "Requirement already satisfied: pytz>=2020.1 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from pandas>=1.2->seaborn->-r requirements.txt (line 1)) (2025.2)\n", "Requirement already satisfied: tzdata>=2022.7 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from pandas>=1.2->seaborn->-r requirements.txt (line 1)) (2025.2)\n", "Requirement already satisfied: ptyprocess>=0.5 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from pexpect>4.3->ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (0.7.0)\n", "Requirement already satisfied: certifi in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from pyogrio>=0.7.2->geopandas->-r requirements.txt (line 2)) (2025.8.3)\n", "Requirement already satisfied: six>=1.5 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from python-dateutil>=2.7->matplotlib!=3.6.1,>=3.4->seaborn->-r requirements.txt (line 1)) (1.17.0)\n", "Requirement already satisfied: executing>=1.2.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from stack_data->ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (2.2.0)\n", "Requirement already satisfied: asttokens>=2.1.0 in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from stack_data->ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (3.0.0)\n", "Requirement already satisfied: pure-eval in /home/gep/tmp/ws3/.venv/lib/python3.12/site-packages (from stack_data->ipython>=6.1.0->ipywidgets->-r requirements.txt (line 3)) (0.2.3)\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "%pip install -r requirements.txt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import required packages." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import ws3\n", "import ws3.forest\n", "import functools\n", "import pandas as pd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set up ``ForestModel`` instance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set some reasonable model parameters." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "base_year = 2020 # base year for the problem\n", "horizon = 10 # number of periods in the simulation horizon\n", "period_length = 10 # period length (in years)\n", "max_age = 1000 # maximum age of a stand (in years)\n", "tvy_name = \"totvol\" # name for total volume yield component\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a new `ForestModel` object and import a model dataset from the `data` directory." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "fm = ws3.forest.ForestModel(model_name=\"tsa24_clipped\",\n", " model_path=\"data/woodstock_model_files_tsa24_clipped\",\n", " base_year=base_year,\n", " horizon=horizon,\n", " period_length=period_length,\n", " max_age=max_age)\n", "fm.import_landscape_section()\n", "fm.import_areas_section(convert_periods_to_years=period_length)\n", "fm.import_yields_section(convert_periods_to_years=period_length)\n", "fm.import_actions_section(convert_periods_to_years=period_length)\n", "fm.import_transitions_section(convert_periods_to_years=period_length)\n", "fm.initialize_areas()\n", "fm.add_null_action()\n", "fm.reset_actions()\n", "\n", "fm.actions[\"harvest\"].is_harvest = True # set harvest action to be a harvest action (needed for `cmp_c_z` to work correctly)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Explore optimization modelling functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define a base scenario (maximizing harvest volume)\n", "\n", "We will start by defining a basic optimization problem, where we simply maximize the sum of harvest volume across all time periods. \n", "\n", "The optimization functions in `ws3` are very flexible, but this flexibility comes at a cost: the user has to define various functions that generate the linear programming (LP) optimization problem matrix coefficients.\n", "\n", "We will start by defining a function `cmp_c_z` (you can use any name you want for coefficient functions) that returns the objective function coefficient for a single column in the matrix. Note that `ws3` is using a Model I formulation of the LP problem, where each column of the matrix represents a unique feasible combination of actions across all time periods (for a given combination of development type and initial age class). Internally `ws3` handles iterating through each development type/age class combination, and generating dynamic programming state trees, where each unique _path_ from the root node of a tree through to a leaf node is represented by a unique column in the matrix. Each path has one node per time period in the simulation horizon.\n", "\n", "The `cmp_c_z` function will be used by `ws3` to generate the coefficients for each column in the matrix. The `cmp_c_z` function takes three arguments: the first argument is a `ForestModel instance`, the second argument is a `path` (a list of node objects), and the third argument is an `expr` (a string expression that evaluates to a number, which we pass to the `ForestModel.compile_product` method).\n", "\n", "See the [ws3 documentation](https://ws3.readthedocs.org) for more information about what types of expressions `ws3` can resolve in this context. For our example we will use a simple expression that multiplies total merchantable volume in the `totvol` yield component by 0.85 to simulate an 85% fibre utilization rate (i.e., 15% of total merchnantable volume is left behind as harvest residues in slash piles or spread across cut blocks).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Model I LP optimization problem can be formulated as follows:\n", "\n", "$$\n", "\\begin{align}\n", " \\text{max} \\quad & \\sum_{i\\in I}\\sum_{j \\in J_{i}}c_{ij}x_{ij} & \\\\\n", " \\text{s.t.} \\quad & \\nonumber \\\\\n", " &(1-\\varepsilon_{p})y_p \\leq \\sum_{i\\in I}\\sum_{j\\in J_{i}}\\mu_{ijpt}x_{ij} \\leq (1+\\varepsilon_{p})y_p, & \\forall p\\in O', t\\in T'_p \\\\\n", " &v^{-}_{ot}\\leq\\sum_{i\\in Z}\\sum_{j\\in J_{i}}\\mu_{ijot}x_{ij} \\leq v^{+}_{ot}, & \\forall o\\in O, t\\in T \\\\\n", " &\\sum_{j\\in J_{i}}x_{ij} = 1, \\forall i\\in Z \\\\\n", " &0 \\leq x_{ij} \\leq 1, & \\forall i \\in Z, j \\in J_{i} \n", "\\end{align}\n", "$$\n", "\n", "where\n", "\n", "$$\n", "\\begin{align*}\n", " I := & \\,\\, \\text{set of spatial zones}\\\\\n", " J_{i} := & \\,\\, \\text{set of available prescriptions for zone $i \\in I$}\\\\\n", " O := & \\,\\, \\text{set of forest outputs}\\\\\n", " O' \\subseteq O := & \\,\\, \\text{set of targeted forest outputs}\\\\\n", " T := & \\,\\, \\text{set of time periods in the planning horizon}\\\\\n", " T'_p \\subseteq O := & \\,\\, \\text{subset of $T$ on which even-flow constraints for output $p \\in O^{\\prime}$ are applied}\\\\\n", " \\varepsilon_{p} := & \\,\\, \\text{admissible level of variation on yield of targeted output $p \\in O^{\\prime}$} \\\\\n", " \\mu_{ijot} := & \\,\\, \\text{quantity of output $o \\in O$ produced in period $t \\in T$ by prescription $j \\in J_{i}$}\\\\\n", " & \\,\\, \\text{in zone $i \\in I$}\\\\\n", " \\mu_{ijpt} := & \\,\\, \\text{quantity of output $p \\in O'$ produced in period $t \\in T$ by prescription $j \\in J_{i}$}\\\\\n", " & \\,\\, \\text{in zone $i \\in I$}\\\\\n", " v^{-}_{ot} := & \\,\\, \\text{lower bound on yield of output $o \\in O$ in period $t \\in T$} \\\\\n", " v^{+}_{ot} := & \\,\\, \\text{upper bound on yield of output $o \\in O$ in period $t \\in T$} \\\\\n", " c_{ij} := & \\,\\, \\text{objective function contribution of prescription $j \\in J_{i}$ in zone $i \\in I$} \\\\\n", " x_{ij} := & \\,\\, \\text{proportion of zone $i \\in Z$ on which prescription $j \\in J_{i}$ is applied} \\\\\n", " y_p = & \\sum_{i\\in I}\\sum_{j\\in J_i}\\mu_{ijpt^R_p}x_{ij}\\text{, i.e. total yield of targeted output $p \\in O'$ at reference period $t^R_p$} \n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "The objective function (1) maximizes the sum of $c_{ij}x_{ij}$ products, which represent yield of a user-defined output---the output can be anything, but will typically be harvest volume or harvested area.\n", "\n", "The variables $x_{ij}$ are linear, with domain $\\{x_{ij} \\in \\mathbb{R}|0 \\leq x_{ij} \\leq 1\\}$ specified in (5).\n", "\n", "The set $O^{\\prime} \\subseteq O$ represents targeted outputs, for which we enforce even-flow constraints. Even-flow constraints are defined in (2), and are expressed in terms of $y_p$, which represents total yield of targeted output $p \\in O^{\\prime}$ in reference period $t^R_p$. Note that even-flow constraints are defined over time periods $T'_p$, which is any subset of $T$ (need not be contiguous)---$T'_p$ can be unique for each output $p \\in O^{\\prime}$.\n", "\n", "General constraints (3) set upper and lower bounds on periodic yield of any output $o \\in O$---we can use these constraints to set minimum or maximum periodic levels of any indicators that we can formulate from data in the model (e.g., set an upper bount on harvest volume in the first period, or set a lower bound on growing stock in the final period).\n", "\n", "Coverage constraints (4) are accounting constraints that force prescriptions to cover the entire zone---doing nothing for the entire planning horizon is considered a prescription that could generate some outputs. This maintains the structure of the model, and these are automatically defined and set by `ws3` (i.e., the user does not need to define these constraints as part of optimization model input)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def cmp_c_z(fm, path, expr):\n", " \"\"\"\n", " Compile objective function coefficient (given ForestModel instance, \n", " root-to-leaf state tree node path, and an expression that resolves to the correct product unit).\n", " \n", " This function iterates over the specified node path in a state tree, \n", " checks if each node corresponds to a harvest operation, and adds up the compiled product values for these nodes.\n", " \n", " Parameters:\n", " fm (ForestModel): The ForestModel instance containing the necessary information.\n", " path (list of Node objects): A list of root-to-leaf nodes forming a path through a state tree.\n", " expr (str): An expression resolving to the correct product unit.\n", " \n", " Returns:\n", " float: The compiled objective function coefficient.\n", " \"\"\"\n", " result = 0. # Initialize the result variable\n", " for t, n in enumerate(path, start=1): # Iterate over the node path with index and value\n", " d = n.data() # Get the data from the current node\n", " if fm.is_harvest(d[\"acode\"]): # Check if the current node corresponds to a harvest operation\n", " result += fm.compile_product(t, expr, d[\"acode\"], [d[\"dtk\"]], d[\"age\"], coeff=False) # Add up the compiled product values\n", " return result # Return the final compiled objective function coefficient\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that a decision variables in the Model I problem represents a proportion of the total area in a given tree that is assigned to a given path (sequence of actions), and that the sum of all decision variables corresponding to a given tree is automatically constrained to sum to 1. As such, the objective function coefficients should represent absolute yield (e.g., $m^3$ for harvest volume, not $m^3 ha^{-1}$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we define the `coeff_funcs` dictionary, which maps matrix row names to the functions that we defined that generate the correct coefficients for each column of the matrix in that row. In this simple we only map row `z` names to the `cmp_c_z` function that we defined above. This tells the `ForestModel.add_problem` method to use the `cmp_c_z` function to compile the objective function coefficients (the objective function is named `z` by default in the `ForestModel.add_problem` method arg signature definition, so if we map a function to name `z` in `coeff_funcs` it will automatically pick that up as the objective function coefficient function).\n", "\n", "Note that the `add_problem` method assumes that all coefficient functions passed to `coeff_funcs` will have exactly two arguments (a `ForestModel` instance and a `Path` instance). If this is not the case then you will need to define a wrapper function that rewrites the function signature (see the example below using `functools.partial`)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "expr = \"0.85 * totvol\"\n", "coeff_funcs = {\"z\":functools.partial(cmp_c_z, expr=expr)}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to call the `ForestModel.add_problem` method, which will build the optimization problem (including the objective function coefficients from the function we provided) and add the problem to the model under the name `'scenario-1'`. The `add_problem` method returns an instance of the `Problem` class that we can use to solve the LP problem and access the optimal solution (which we will compile into a `ws3` action schedule so we can simulate it out and plot results in a human-readable way). \n", "\n", "We pass `None` value to the `cflw_e` and `cgen_data` arguments because there are no flow or general constraints in this problem.\n", "\n", "We pass `'null'` and `'harvest'` values to the `acodes` argument because these are the action we want to include in the generation of the Model 1 state trees in this problem. This allows us to define complex `ForestModel` base models, which potentially define a wide range of actions, without being required to always include all the defined actions in all the optimization problems. It is required to define a _null_ action in the `ForestModel` instance before calling the `add_problem` method (the `ForestModel.add_null_action` method makes it easy to correctly define a null action), and you must include the name of the null action in the list of action codes passed to the `acodes` argument of the `add_problem` method. Including a null action allows the optimization problem to simulate doing nothing (i.e., just letting the forest grow). If you want to include all the actions defined in the `ForestModel`, pass `None` value to the `acodes` argument of the `add_problem` method.\n", "\n", "We pass `None` to the `mask` argument because we want to include all the development types in this problem. You can filter which developments are included in a given problem by passing a valid mask (either in tuple format or as a Woodstock-style mask string). See the [ws3 documentation](https://ws3.readthedocs.org) for more information about themes and masks (including how to define aggregate landscape codes to match multiple values within one theme in a mask).\n", "\n", "Note also that the objective function sense gets defined at this point by passing the `sense` argument to the `add_problem` method. The `sense` argument can be either `ws3.opt.SENSE_MINIMIZE` or `ws3.opt.SENSE_MAXIMIZE`." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "problem = fm.add_problem(name=\"scenario-1\",\n", " coeff_funcs=coeff_funcs, \n", " cflw_e=None, \n", " cgen_data=None, \n", " acodes=(\"null\", \"harvest\"),\n", " sense=ws3.opt.SENSE_MAXIMIZE, \n", " mask=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have a `Problem` object instance, we can call the `solve` method to try to solve the problem. \n", "\n", "This will use the default bindings to the open source HiGHS solver (see the [HiGHS documentation](https://www.highs.dev/) for details on HiGHS) to find a solution that maximizes the objective. " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running HiGHS 1.11.0 (git hash: 364c83a): Copyright (c) 2025 HiGHS under MIT licence terms\n", "LP has 32 rows; 305 cols; 305 nonzeros\n", "Coefficient ranges:\n", " Matrix [1e+00, 1e+00]\n", " Cost [2e+01, 4e+04]\n", " Bound [1e+00, 1e+00]\n", " RHS [1e+00, 1e+00]\n", "Presolving model\n", "24 rows, 115 cols, 115 nonzeros 0s\n", "0 rows, 0 cols, 0 nonzeros 0s\n", "Presolve : Reductions: rows 0(-32); columns 0(-305); elements 0(-305) - Reduced to empty\n", "Solving the original LP from the solution after postsolve\n", "Using EKK primal simplex solver\n", " Iteration Objective Infeasibilities num(sum)\n", " 0 -2.2455249290e+05 Pr: 0(0); Du: 222(1.63573e+06) 0s\n", " 19 -2.2455284502e+05 Pr: 0(0) 0s\n", "Required 19 simplex iterations after postsolve\n", "Model status : Optimal\n", "Simplex iterations: 19\n", "Objective value : -2.2455284502e+05\n", "P-D objective error : 6.4803823015e-17\n", "HiGHS run time : 0.00\n" ] } ], "source": [ "problem.solve(verbose=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use the `Problem.status` method to check the status of the solution (i.e., to confirm that the problem was solved to optimality)." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Hooray! We have an optimal solution!\n" ] } ], "source": [ "if (problem.status() == ws3.opt.STATUS_OPTIMAL):\n", " print(\"Hooray! We have an optimal solution!\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use the `Problem.z` method to check the value of the objective function at the optimal solution." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The optimal value of the objective function is: 224552.84502455403\n" ] } ], "source": [ "print(\"The optimal value of the objective function is: {}\".format(problem.z()))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This matches the value reported by HiGHS in the solver `stdout` console output in a previous step." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we need to extract the optimal sequence of actions from the `Problem` instance and compile it into a format that can be used as input to the `ForestModel.apply_schedule` method. We can use the `ForestModel.compile_schedule` method to do this. Internally this dispatches the call to private method `ForestModel._cmp_sch_m1`, which is aware of the Model I optimization problem structure and correctly walks through all the `Node` objects in all the `Path` objects in all the `Tree` objects in the `Problem` instance that was passed in as an argument, extracting the required data (i.e., `DevelopmentType` key, age, actionned area, action code, period, existing/future `DevelopmentType` status) from the `data` dictionary of each `Node` object and assembling this information into a correctly formatted list of tuples that can be consumed by `ForestModel.apply_action`." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 85,\n", " 282.1296355046733,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 91,\n", " 73.1021561503533,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 93,\n", " 28.37956666951611,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 95,\n", " 94.94675966211176,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 105,\n", " 32.175418531537815,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 113,\n", " 4.184826329641321,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 115,\n", " 50.030816858894816,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 125,\n", " 78.16612132001225,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 135,\n", " 72.24421919373785,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 145,\n", " 96.38442685503642,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 153,\n", " 9.591469412607397,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 155,\n", " 34.32629241113743,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402000', '100', '2402000'),\n", " 165,\n", " 0.638005468748551,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2402002'),\n", " 93,\n", " 48.21816452980633,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2402002'),\n", " 95,\n", " 33.89498244313859,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2402002'),\n", " 115,\n", " 3.195378954654358,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2403000', '100', '2403000'),\n", " 93,\n", " 14.811643286926836,\n", " 'harvest',\n", " 1,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2402002'),\n", " 88,\n", " 32.64168183055375,\n", " 'harvest',\n", " 2,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2403002', '204', '2403002'),\n", " 83,\n", " 2.243990590272984,\n", " 'harvest',\n", " 2,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 168,\n", " 103.76740323520823,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2401002'),\n", " 170,\n", " 4.173147018452507,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 282.1296355046733,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 73.1021561503533,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 28.37956666951611,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 94.94675966211176,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 32.175418531537815,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 4.184826329641321,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 50.030816858894816,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 78.16612132001225,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 72.24421919373785,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 96.38442685503642,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 9.591469412607397,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 90,\n", " 34.32629241113743,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2401002', '204', '2421002'),\n", " 110,\n", " 0.422054121206099,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402000', '100', '2422000'),\n", " 90,\n", " 0.638005468748551,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2422002'),\n", " 80,\n", " 32.64168183055375,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2422002'),\n", " 90,\n", " 48.21816452980633,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2422002'),\n", " 90,\n", " 33.89498244313859,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2402002', '204', '2422002'),\n", " 90,\n", " 3.195378954654358,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2403000', '100', '2423000'),\n", " 90,\n", " 14.811643286926836,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2403002', '204', '2423002'),\n", " 80,\n", " 2.243990590272984,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2403002', '204', '2423002'),\n", " 99,\n", " 59.81429119367274,\n", " 'harvest',\n", " 10,\n", " '_existing'),\n", " (('tsa24_clipped', '1', '2403002', '204', '2423002'),\n", " 108,\n", " 32.366198551219505,\n", " 'harvest',\n", " 10,\n", " '_existing')]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "schedule = fm.compile_schedule(problem)\n", "schedule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the schedule is just a list of 6-element tuples. \n", "\n", "The first element in each tuple is the `DevelopmentType` key (a tuple of landscape theme strings).\n", "\n", "The other elements in the tuple specify the age, actionned area, action code, period, and existing/future `DevelopmentType` status. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we that we have a compiled schedule, we can use it as input for the `ForestModel.apply_schedule` method using the following parameters:\n", "- `force_integral_area`: Set to `False` to avoid forcing integral area calculations\n", "- `override_operability`: Set to `False` to avoid overriding operability calculations\n", "- `fuzzy_age`: Set to `False` to disable fuzzy age calculations\n", "- `recourse_enabled`: Set to `False` to disable recourse-enabled mode\n", "- `verbose`: Set to `False` for quiet operation\n", "- `compile_c_ycomps`: Set to `True` to enable compilation of complex yield components\n", "\n", "See the [ws3 documentation](https://ws3.readthedocs.org) for the `ForestModel.apply_schedule` and `ForestModel.apply_action` methods for more information about these and other parameters that are available to controls exactly how actions get applied to the model. Ultimately, the `apply_schedule` method is a wrapper around the `apply_action` method. The `apply_action` method will handle the details of applying individual actions to the model. One of the main functions of `ws3` (and any forest estate model) is to simulate forest response to sequences of actions interleaved with periods of growth. Given the core role of the `apply_action` method, it is not surprising that it has a large number of parameters that can be used to control exactly how actions are applied.\n", "\n", "This effectively resets the current schedule and model state to the initial state, and runs a sequential multi-period simulation of alternating actions and periods of growth. " ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "fm.apply_schedule(\n", " schedule, \n", " force_integral_area=False, \n", " override_operability=False,\n", " fuzzy_age=False,\n", " recourse_enabled=False,\n", " verbose=False,\n", " compile_c_ycomps=True\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point, our `ForestModel` object has simulated sequential application of the optimal schedule of activities we created in a previous step." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we compile some relevant performance indicators to summarize what is going in this optimal solution (i.e., periodic harvested area, periodic harvest volume, periodic growing stock).\n", "\n", "First, we compile the harvested area and volume indicators. We use the `ForestModel.ompile_product` method to compile these indicators. The first argument is the period for which we want to compile a product output, while the second argument is the expression that will be resolved by `ws3` to compute a yield coefficient. The `compile_product` method compiles the product of actionned area and whatever yield coefficient the provided expression resolves to (expressed in _units of output product per unit of area_, e.g., $m^3/ha$ for harvested volume or $ha/ha$ for harvested area). So, when we use an expression value of `1.` for the second argument of the `compile_product` method, we are telling `ws3` that we want to _pass through_ the actionned area indicator that it is calculating internally. We use the `acode` argument to specify which action code to use (i.e., `harvest` to compile harvested area and harvested volume indicators).\n", "\n", "We use the `ForestModel.inventory` method to compile the growing stock indicator. The first argument is the period for which we want to compile an inventory output, while the second argument is the expression that will be resolved by `ws3` to compute a yield coefficient (in this case we use the total volume of trees multiplied by an assumed utilization rate of 85%). \n", "\n", "See the [ws3 documentation](https://ws3.readthedocs.org) for more information on expression evaluation, and what types of expressions are supported by `ws3` that can be provided as input to `compile_product` and `inventory` method calls.\n", "\n", "Once we have compiled our raw indictor outputs into lists (with one value period, including a data column of period index values), we stuff all of this into a dictionary and create a new `pandas.DataFrame` to hold our compiled data. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
| \n", " | period | \n", "oha | \n", "ohv | \n", "ogs | \n", "
|---|---|---|---|---|
| 0 | \n", "1 | \n", "956.419884 | \n", "102585.020390 | \n", "30405.125309 | \n", "
| 1 | \n", "2 | \n", "34.885672 | \n", "4798.494437 | \n", "28607.304394 | \n", "
| 2 | \n", "3 | \n", "0.000000 | \n", "0.000000 | \n", "32154.392350 | \n", "
| 3 | \n", "4 | \n", "0.000000 | \n", "0.000000 | \n", "38166.885368 | \n", "
| 4 | \n", "5 | \n", "0.000000 | \n", "0.000000 | \n", "47379.394949 | \n", "
| 5 | \n", "6 | \n", "0.000000 | \n", "0.000000 | \n", "59578.656570 | \n", "
| 6 | \n", "7 | \n", "0.000000 | \n", "0.000000 | \n", "77531.535059 | \n", "
| 7 | \n", "8 | \n", "0.000000 | \n", "0.000000 | \n", "103522.146540 | \n", "
| 8 | \n", "9 | \n", "0.000000 | \n", "0.000000 | \n", "134548.020901 | \n", "
| 9 | \n", "10 | \n", "1191.848650 | \n", "117169.330197 | \n", "25643.073440 | \n", "