# Alkylation Process Optimization (PROCESS)#

`process.py`

```"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_process.html
## 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
from gamspy import Equation
from gamspy import Model
from gamspy import Problem
from gamspy import Sense
from gamspy import Variable
from gamspy.math import sqr

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

# 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("Profit in model 'process': {:.2f}".format(profit.records.level[0]))
rproc.solve()
print("Profit in model 'rproc': {:.2f}".format(profit.records.level[0]))

if __name__ == "__main__":
main()
```