Alkylation Process Optimization (PROCESS)

process.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_process.html
## LICENSETYPE: Demo
## MODELTYPE: NLP
## KEYWORDS: nonlinear programming, alkylation process, chemical engineering


Alkylation Process Optimization (PROCESS)

Optimization of a alkylation process.


Bracken, J, and McCormick, G P, Chapter 4. In Selected Applications
of Nonlinear Programming. John Wiley and Sons, New York, 1968.
"""

from __future__ import annotations

import os

from gamspy import Container, Equation, Model, Problem, Sense, Variable
from gamspy.math import sqr


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
    )

    # Variables
    olefin = Variable(
        m, name="olefin", type="positive", description="olefin feed (bpd)"
    )
    isor = Variable(
        m, name="isor", type="positive", description="isobutane recycle (bpd)"
    )
    acid = Variable(
        m,
        name="acid",
        type="positive",
        description="acid addition rate (1000lb per day)",
    )
    alkylate = Variable(
        m, name="alkylate", type="positive", description="alkylate yield (bpd)"
    )
    isom = Variable(
        m, name="isom", type="positive", description="isobutane makeup (bpd)"
    )
    strength = Variable(
        m,
        name="strength",
        type="positive",
        description="acid strength (weight pct)",
    )
    octane = Variable(
        m, name="octane", type="positive", description="motor octane number"
    )
    ratio = Variable(
        m,
        name="ratio",
        type="positive",
        description="isobutane makeup to olefin ratio",
    )
    dilute = Variable(
        m, name="dilute", type="positive", description="acid dilution factor"
    )
    f4 = Variable(
        m, name="f4", type="positive", description="f-4 performance number"
    )

    profit = Variable(m, name="profit")
    rangey = Variable(m, name="rangey")
    rangem = Variable(m, name="rangem")
    ranged = Variable(m, name="ranged")
    rangef = Variable(m, name="rangef")

    # Equations
    yield1 = Equation(
        m, name="yield1", description="alkylate yield definition"
    )
    rngyield = Equation(
        m, name="rngyield", description="ranged alkylate yield definition"
    )
    makeup = Equation(
        m, name="makeup", description="isobutane makeup definition"
    )
    sdef = Equation(m, name="sdef", description="acid strength definition")
    motor = Equation(m, name="motor", description="motor octane number")
    rngmotor = Equation(
        m, name="rngmotor", description="ranged motor octane number"
    )
    drat = Equation(m, name="drat", description="isobutane to olefin ratio")
    ddil = Equation(m, name="ddil", description="dilution definition")
    rngddil = Equation(
        m, name="rngddil", description="ranged dilution definition"
    )
    df4 = Equation(m, name="df4", description="f-4 definition")
    rngdf4 = Equation(m, name="rngdf4", description="ranged f-4 definition")
    dprofit = Equation(m, name="dprofit", description="profit definition")

    yield1[...] = alkylate == olefin * (
        1.12 + 0.13167 * ratio - 0.00667 * sqr(ratio)
    )
    makeup[...] = alkylate == olefin + isom - 0.22 * alkylate
    sdef[...] = acid == alkylate * dilute * strength / (98 - strength) / 1000
    motor[...] = octane == 86.35 + 1.098 * ratio - 0.038 * sqr(
        ratio
    ) - 0.325 * (89 - strength)
    drat[...] = ratio == (isor + isom) / olefin
    ddil[...] = dilute == 35.82 - 0.222 * f4
    df4[...] = f4 == -133 + 3 * octane
    dprofit[...] = (
        profit
        == 0.063 * alkylate * octane
        - 5.04 * olefin
        - 0.035 * isor
        - 10 * acid
        - 3.36 * isom
    )
    rngyield[...] = rangey * alkylate == olefin * (
        1.12 + 0.13167 * ratio - 0.00667 * sqr(ratio)
    )
    rngmotor[...] = rangem * octane == 86.35 + 1.098 * ratio - 0.038 * sqr(
        ratio
    ) - 0.325 * (89 - strength)
    rngddil[...] = ranged * dilute == 35.82 - 0.222 * f4
    rngdf4[...] = rangef * f4 == -133 + 3 * octane

    # Define Models
    process = Model(
        m,
        name="process",
        equations=[yield1, makeup, sdef, motor, drat, ddil, df4, dprofit],
        problem=Problem.NLP,
        sense=Sense.MAX,
        objective=profit,
    )
    rproc = Model(
        m,
        name="rproc",
        equations=[
            rngyield,
            makeup,
            sdef,
            rngmotor,
            drat,
            rngddil,
            rngdf4,
            dprofit,
        ],
        problem=Problem.NLP,
        sense=Sense.MAX,
        objective=profit,
    )

    rangey.lo[...] = 0.9
    rangey.up[...] = 1.1
    rangey.l[...] = 1
    rangem.lo[...] = 0.9
    rangem.up[...] = 1.1
    rangem.l[...] = 1
    ranged.lo[...] = 0.9
    ranged.up[...] = 1.1
    ranged.l[...] = 1
    rangef.lo[...] = 0.9
    rangef.up[...] = 1.1
    rangef.l[...] = 1

    strength.lo[...] = 85
    strength.up[...] = 93
    octane.lo[...] = 90
    octane.up[...] = 95
    ratio.lo[...] = 3
    ratio.up[...] = 12
    dilute.lo[...] = 1.2
    dilute.up[...] = 4
    f4.lo[...] = 145
    f4.up[...] = 162
    olefin.lo[...] = 10
    olefin.up[...] = 2000
    isor.up[...] = 16000
    acid.up[...] = 120
    alkylate.up[...] = 5000
    isom.up[...] = 2000

    olefin.l[...] = 1745
    isor.l[...] = 12000
    acid.l[...] = 110
    alkylate.l[...] = 3048
    isom.l[...] = 1974
    strength.l[...] = 89.2
    octane.l[...] = 92.8
    ratio.l[...] = 8
    dilute.l[...] = 3.6
    f4.l[...] = 145
    profit.l[...] = 872

    process.solve()
    print(f"Profit in model 'process': {profit.records.level[0]:.2f}")
    rproc.solve()
    print(f"Profit in model 'rproc': {profit.records.level[0]:.2f}")


if __name__ == "__main__":
    main()