Metadata-Version: 2.4
Name: pypoly-reactor
Version: 0.1.3
Summary: Polymerization reactor modeling package using the method of moments
Author: Dabin Yang
License-Expression: MIT
Keywords: polymerization,reactor modeling,method of moments,chemical engineering
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Provides-Extra: plot
Requires-Dist: matplotlib>=3.7; extra == "plot"
Dynamic: license-file

# PyPoly Reactor + RL

![PyPolymer logo](src/pypoly/assets/pypolymer-logo.png)

`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:

```text
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

```bash
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:

```bash
pip install .
```

Then import:

```python
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:

```python
import pypoly

logo = pypoly.logo_path()
print(logo)
```

---

# Part A — Reactor simulation (`PolyRxn`)

## Batch / CSTR

```python
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}`):

```text
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.

```python
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):

```text
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

```python
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

```python
from pypoly import run_cstr_rl

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

## LDPE PFR RL

```python
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:

```text
  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):

```python
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:

```python
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:

  ```text
  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

```bash
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.
