Antalya Forestry Model - Steady State (TFORSS)#

tforss.py tforss.gdx

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_tforss.html
## LICENSETYPE: Demo
## MODELTYPE: LP
## DATAFILES: tforss.gdx
## KEYWORDS: linear programming, forestry, scenario analysis, investment planning, forest management planning


Antalya Forestry Model - Steady State (TFORSS)

This model finds the best management plan for new forests in a steady state
condition.


Bergendorff, H, Glenshaw, P, and Meeraus, A, The Planning of Investment
Programs in the Paper Industry. Tech. rep., The World Bank, 1980.
"""

from __future__ import annotations

import os
from pathlib import Path

import numpy as np

from gamspy import Container
from gamspy import Equation
from gamspy import Model
from gamspy import Ord
from gamspy import Parameter
from gamspy import Sense
from gamspy import Set
from gamspy import Sum
from gamspy import Variable
from gamspy.math import power
from gamspy.math import Round


def main():
    cont = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
        load_from=str(Path(__file__).parent.absolute()) + "/tforss.gdx",
    )

    # Sets
    c, cf, cl, s, k, at, p, m = cont.getSymbols(
        ["c", "cf", "cl", "s", "k", "at", "p", "m"]
    )

    # Parameters
    scd, land, ymf, a, b, pc, pd, nu, age = cont.getSymbols([
        "scd",
        "land",
        "ymf",
        "a",
        "b",
        "pc",
        "pd",
        "nu",
        "age",
    ])

    # Scalar
    mup, muc, life, rho = cont.getSymbols([
        "mup",
        "muc",
        "life",
        "rho",
    ])

    age[at] = 10 * Ord(at)

    # Model Definition #

    # Equation
    lbal = Equation(
        cont,
        name="lbal",
        domain=cl,
        description="log balances",
    )
    bal = Equation(
        cont,
        name="bal",
        domain=c,
        description="material balances of wood processing",
    )
    cap = Equation(
        cont,
        name="cap",
        domain=m,
        description="wood processing capacities",
    )
    landc = Equation(
        cont,
        name="landc",
        domain=[s, k],
        description="land availability constraint",
    )
    ainvc = Equation(cont, name="ainvc", description="investment cost")
    aproc = Equation(cont, name="aproc", description="process cost")
    asales = Equation(cont, name="asales", description="sales revenue")
    acutc = Equation(cont, name="acutc", description="cutting cost")
    aplnt = Equation(cont, name="aplnt", description="planting cost")

    # Variable
    v = Variable(
        cont,
        name="v",
        type="positive",
        domain=[s, k, at],
        description="management of new forest   (1000ha per year)",
    )
    r = Variable(
        cont,
        name="r",
        domain=c,
        description="supply of logs to industry (1000m3 per year)",
    )
    z = Variable(
        cont,
        name="z",
        type="positive",
        domain=p,
        description="process level        (1000m3 input per year)",
    )
    h = Variable(
        cont,
        name="h",
        domain=m,
        description="capacity             (1000m3 input per year)",
    )
    x = Variable(
        cont,
        name="x",
        type="positive",
        domain=c,
        description="final shipments        (1000 units per year)",
    )
    phik = Variable(
        cont,
        name="phik",
        description="investment cost           (1000us$ per year)",
    )
    phir = Variable(
        cont,
        name="phir",
        description="process cost              (1000us$ per year)",
    )
    phix = Variable(
        cont,
        name="phix",
        description="sales revenue             (1000us$ per year)",
    )
    phil = Variable(
        cont,
        name="phil",
        description="cutting cost              (1000us$ per year)",
    )
    phip = Variable(
        cont,
        name="phip",
        description="planting cost             (1000us$ per year)",
    )

    lbal[cl] = r[cl] == Sum([s, k, at], ymf[at, k, s, cl] * v[s, k, at])

    bal[c] = Sum(p, a[c, p] * z[p]) + r[c].where[cl[c]] >= x[c].where[cf[c]]

    cap[m] = Sum(p, b[m, p] * z[p]) == h[m]

    landc[s, k] = Sum(at, v[s, k, at] * age[at]) <= land[s] * scd[k]

    ainvc[...] = phik == rho / (1 - power((1 + rho), (-life))) * Sum(
        m, nu[m] * h[m]
    )

    aproc[...] = phir == Sum(p, pc[p] * z[p])

    asales[...] = phix == Sum(cf, pd[cf] * x[cf])

    acutc[...] = phil == muc * Sum(cl, r[cl])

    aplnt[...] = phip == mup * Sum(
        [s, k, at], v[s, k, at] * (1 + rho) ** age[at]
    )

    # Objective Function; total benefits             (discounted cost)
    benefit = phix - phik - phir - phil - phip

    # Model definition
    forest = Model(
        cont,
        name="forest",
        equations=cont.getEquations(),
        problem="LP",
        sense=Sense.MAX,
        objective=benefit,
    )

    # Case Selection and Report Definitions
    rhoset = Set(
        cont, name="rhoset", records=["rho-03", "rho-05", "rho-07", "rho-10"]
    )

    landcl = Parameter(
        cont, name="landcl", domain=[s, k], description="clean level of landc"
    )
    rep = Parameter(
        cont,
        name="rep",
        domain=[cl, rhoset],
        description="summary report on log supply      (1000m3 per year)",
    )
    reprp = Parameter(
        cont,
        name="reprp",
        domain=[s, k, rhoset],
        description="summary report on rotation period           (years)",
    )
    repsp = Parameter(
        cont,
        name="repsp",
        domain=[s, k, rhoset],
        description="summary report on shadow price of land (us$ per ha)",
    )
    rhoval = Parameter(
        cont,
        name="rhoval",
        domain=rhoset,
        records=np.array([0.03, 0.05, 0.07, 0.1]),
    )

    for case, _ in rhoset.records.itertuples(index=False):
        rho[...] = rhoval[case]
        forest.solve()
        landcl[s, k] = Round(landc.l[s, k], 3)
        rep[cl, case] = r.l[cl]
        reprp[s, k, case] = (landcl[s, k] / Sum(at, v.l[s, k, at])).where[
            landcl[s, k]
        ]
        repsp[s, k, case] = landc.m[s, k]

    print("Summary report on log supply (1000m3 per year):")
    print(rep.pivot().round(3))
    print()

    print("Summary report on rotation period (years)  :")
    print(reprp.pivot().round(3))
    print()

    print("Summary report on shadow price of land (us$ per ha)  :")
    print(repsp.pivot().round(3))


if __name__ == "__main__":
    main()