pypoly-reactor 0.1.3

Polymerization reactor modeling package using the method of moments

PyPoly Reactor + RL

pypoly-reactor 0.1.3 packages a single-file Python module that combines:

1. a polymerization reactor simulator (PolyRxn: Batch / CSTR / tubular PFR), and
2. reinforcement-learning (RL) operating-condition optimisation for all three reactors (Batch, CSTR, and LDPE PFR), using the same REINFORCE method.

The implementation is packaged as:

src/pypoly/pypoly.py

Part	Contents
1	PolyRxn simulator — Batch, CSTR, tubular PFR (method of moments). PFR has optional beta-scission.
2	LDPE PFR RL — the research workflow (run_pfr, compute_reward, GaussianMLP, REINFORCE, exploitation, N=100 validation) via run_pfr_rl().
2B	Batch / CSTR RL — the same REINFORCE method/policy applied to the lumped reactors via run_batch_rl() / run_cstr_rl().

Importing the module never starts training. The module needs only NumPy + SciPy; matplotlib is optional (plot=True).

---

Installation

conda create -n poly python=3.10 numpy scipy
conda activate poly
pip install pypoly-reactor

For local source distribution testing, run this from the project root instead:

pip install .

Then import:

from pypoly import PolyRxn                                # simulation
from pypoly import run_batch_rl, run_cstr_rl, run_pfr_rl  # optimisation

The package also includes the PyPolymer logo:

import pypoly

logo = pypoly.logo_path()
print(logo)

---

Part A — Reactor simulation (PolyRxn)

Batch / CSTR

from pypoly import PolyRxn

model = PolyRxn(kd=0.8, kp=15.0, ktc=0.02, ktd=0.02,
                ktrm=15e-3, ktrp=10e-3, kca=5e-3, f=0.8,
                dH_p=-100 * 1000, Mw_mono=28.05)

R_gas = 8.3145; Mw_mono = 28.05
P_assu = 3000 * 1e5; T_assu = 150 + 273.15
phys = dict(Tc_const=200 + 273.15,
            rho=Mw_mono * (P_assu / R_gas / T_assu),
            Cp_ass=42.9 / Mw_mono, U_heat=400, D=0.05)

# Batch
model.set_batch_params(**phys)
res = model.run_batch(mono_0=8000.0, ini_0=50.0, CTA_0=50.0,
                      T_0=200 + 273.15, t_end=20.0, dt=0.01)
print("Batch:", res["Mn_final"], res["Mw_final"], res["PDI_final"])

# CSTR (adds residence time tau)
model.set_cstr_params(**phys)
res = model.run_cstr(mono_0=8000.0, ini_0=50.0, CTA_0=50.0,
                     T_0=200 + 273.15, tau=5.0, t_end=20.0, dt=0.01)
print("CSTR: ", res["Mn_final"], res["Mw_final"], res["PDI_final"])

The kinetic constants in these examples are illustrative test values chosen so the code runs quickly, not a calibrated LDPE data set.

PFR (with optional beta-scission)

set_pfr_params gained two optional arguments for beta-scission (Pn -> P_{n/2} + D_{n/2}):

A_kbs   beta-scission pre-exponential   (default 0.0 -> disabled)
Ea_kbs  beta-scission activation energy (default 0.0)

With the defaults the PFR behaves exactly as before; passing the research values (A_kbs=1e11, Ea_kbs=130000) activates scission.

import numpy as np
from scipy.integrate import solve_ivp
from pypoly import PolyRxn

model = PolyRxn(kd=0.8, kp=15.0, ktc=0.02, ktd=0.02,
                ktrm=15e-3, ktrp=10e-3, kca=5e-3, f=0.8,
                dH_p=-100 * 1000, Mw_mono=28.05)

model.set_pfr_params(
    L=10.0, D=0.05, D_j=0.08, v=0.5, v_c=0.5,
    rho_pfr=2.4e6, Cp_pfr=42.9 / 28.05, rho_c=1000.0, Cp_c=4180.0,
    A_kd=0.8, Ea_kd=0.0, A_kp=15.0, Ea_kp=0.0, A_ktc=0.02, Ea_ktc=0.0,
    A_ktd=0.02, Ea_ktd=0.0, A_ktrm=15e-3, Ea_ktrm=0.0,
    A_ktrp=10e-3, Ea_ktrp=0.0, A_kca=5e-3, Ea_kca=0.0,
    P_ref_bar=3000.0, DV_kp=0.0,
    # A_kbs=1e11, Ea_kbs=130000.0,   # uncomment to enable beta-scission
)

N, NV = 50, 11
y0 = np.zeros((N, NV))
y0[:, 6] = 50.0; y0[:, 7] = 8000.0; y0[:, 8] = 50.0
y0[:, 9] = 200 + 273.15; y0[:, 10] = 200 + 273.15

ode = model.pfr_odes(ini_0=50.0, mono_0=8000.0, CTA_0=50.0,
                     T_0=200 + 273.15, Tc_in=200 + 273.15,
                     h_r=0.5, h_j=0.5, N=N)
jac = model.build_jac_sparsity(N, NV)
sol = solve_ivp(ode, (0.0, 5.0), y0.ravel(), method="BDF", jac_sparsity=jac)

print("Success:", sol.success)
print("Final state shape:", sol.y[:, -1].reshape(N, NV).shape)

PolyRxn PFR state vector (NV = 11 per node):

0-2 : live-chain moments   lambda0, lambda1, lambda2
3-5 : dead-chain moments   mu0, mu1, mu2
6   : initiator    7 : monomer    8 : CTA
9   : reactor temperature    10 : coolant temperature

---

Part B — Reinforcement-learning optimisation

A Gaussian-MLP policy is trained with REINFORCE (warm-start, moving-average baseline, Adam) to choose operating conditions that hit product targets. Every run does training -> greedy exploitation -> validation. The Batch, CSTR, and PFR runs all use the same method and policy; only the simulator, levers, and targets differ.

Batch RL

from pypoly import run_batch_rl

out = run_batch_rl(n_eps=120)        # optimise [ini_0, CTA_0, t_end]
print(out["best_oc"], out["valid_result"], out["valid_ok"])

CSTR RL

from pypoly import run_cstr_rl

out = run_cstr_rl(n_eps=120)         # optimise [ini_0, CTA_0, tau]

LDPE PFR RL

from pypoly import run_pfr_rl
# Heavier than Batch/CSTR (a stiff ODE solve per step). Start small.
out = run_pfr_rl(n_eps=150)          # optimise [T0, ini0, Tc_in, U_heat, cta0]

Each call prints a validation table. For example run_batch_rl(n_eps=120, seed=0) gives:

  Results vs targets   [validation: high-fidelity validated]
    X   =  42.22 %       [25–60]   OK
    Tpk =  292.3 °C      [200–300]  OK
    Mn  =  21.65 kg/mol  [12–22]   OK
    Mw  =  44.62 kg/mol  [22–50]  OK
    PDI =   2.06          [1.5–2.5]  OK

Optional learning-curve plot (requires matplotlib):

from pypoly import run_batch_rl
out = run_batch_rl(n_eps=120, plot=True)   # saves rl_batch_results.png

Levers and temperature

Reactor	Levers (action)	Dim
Batch	ini_0, CTA_0, t_end	3
CSTR	ini_0, CTA_0, tau	3
LDPE PFR	T0, ini_0, Tc_in, U_heat, cta0	5

• CTA is the primary Mn lever (more CTA -> lower Mn).
• In Batch/CSTR the rate constants are temperature-independent (only the PFR uses Arrhenius), so feed temperature is not an action. It is still scored: the exotherm raises the temperature with conversion and the reward applies a runaway-safety penalty, so the agent trades conversion against temperature.
• In the PFR, T0 and U_heat are genuine levers.

Configuring targets

Targets are dictionaries you can edit before training. Each entry is (low, high), except temperature which is (low, high, runaway) in degrees C:

import pypoly
pypoly.BATCH_TARGETS = {
    "X":   (0.30, 0.55),
    "T":   (200.0, 290.0, 320.0),
    "Mn":  (15.0, 20.0),    # kg/mol
    "Mw":  (25.0, 45.0),    # kg/mol
    "PDI": (1.8, 2.4),
}
out = pypoly.run_batch_rl(n_eps=120)

---

Model assumptions (please read)

• Empirical high-temperature suppression term. Both PFR formulations include a term of the form RT -= max(0, T - T_safe) * k (T_safe = 573.15 K in the plain PolyRxn PFR, 593.15 K in the RL PFR). This is retained from the original research code for numerical/operational safety and is not a first-principles heat-transfer mechanism — do not read the PFR peak temperature as being limited by jacket heat transfer alone.

• Batch/CSTR rate constants are temperature-independent. kd, kp, ktc, ktd, ktrm, ktrp, kca are constants in the Batch/CSTR models; only the PFR uses Arrhenius kinetics. So in Batch/CSTR the temperature evolves but does not feed back into the reaction rates.

• Two independent PFR implementations with different state ordering. The general PolyRxn PFR and the RL PFR (run_pfr) are separate; their per-node state vectors are ordered differently:

  PolyRxn.pfr_odes : [l0,l1,l2, m0,m1,m2, ini, mono, CTA, T,   Tc ]   (CTA=8, T=9)
  RL run_pfr       : [l0,l1,l2, m0,m1,m2, ini, mono, T,   Tc,  CTA]   (T=8, CTA=10)
  

  Keep this in mind if you index the raw solver output directly.

---

Command-line demo

python -m pypoly              # Batch + CSTR optimisation (fast)
python -m pypoly pfr          # LDPE PFR optimisation (slower, minutes)
python -m pypoly all          # all three

---

Notes

• If using the installed package, import from pypoly. If using the raw single-file module, keep pypoly.py in the same directory as your script or on PYTHONPATH.
• PolyRxn temperatures are in Kelvin; RL tables report peak temperature in degrees C. Mn/Mw are g/mol from PolyRxn, kg/mol in RL output.
• Beta-scission is off by default in PolyRxn (A_kbs=0), so existing simulations are unchanged. The PFR RL workflow uses its own simulator, which already includes beta-scission with the research parameters.
• Reproducibility. Every RL entry point takes a seed argument (run_batch_rl(seed=0), run_cstr_rl(seed=0), run_pfr_rl(seed=0)); the same seed and episode count reproduce the same result.
• Importing the module never starts training; call the run_*_rl() functions explicitly. Batch/CSTR finish in seconds; the PFR run takes minutes.

---

License

Use this file according to the license or distribution agreement provided by the project owner.