"""Pure Python backend for Tennessee Eastman Process simulation. This module provides a complete Python implementation of the TEP simulator, replicating the original Fortran code (teprob.f) without any Fortran dependency. The implementation faithfully reproduces: - All common blocks as Python data structures - The TEINIT initialization routine - The TEFUNC derivative function - All helper subroutines (TESUB1-8) - The linear congruential random number generator - Disturbance random walks and measurement noise This enables running TEP simulations in environments where Fortran compilation is not available, while maintaining statistical similarity to the original. """ import numpy as np from dataclasses import dataclass, field from typing import Optional, Tuple import copy # ============================================================================= # Common Block Data Structures # ============================================================================= @dataclass class ConstBlock: """Physical constants common block (/CONST/). Contains thermodynamic properties for the 8 chemical components: - Antoine coefficients for vapor pressure - Enthalpy coefficients (liquid and gas) - Density coefficients - Heat of vaporization - Molecular weights """ # Molecular weights (kg/kmol) xmw: np.ndarray = field(default_factory=lambda: np.array([ 2.0, 25.4, 28.0, 32.0, 46.0, 48.0, 62.0, 76.0 ])) # Antoine vapor pressure: ln(P) = AVP + BVP/(T + CVP) # Components 1-3 (A,B,C) are non-condensable gases (coefficients = 0) avp: np.ndarray = field(default_factory=lambda: np.array([ 0.0, 0.0, 0.0, 15.92, 16.35, 16.35, 16.43, 17.21 ])) bvp: np.ndarray = field(default_factory=lambda: np.array([ 0.0, 0.0, 0.0, -1444.0, -2114.0, -2114.0, -2748.0, -3318.0 ])) cvp: np.ndarray = field(default_factory=lambda: np.array([ 0.0, 0.0, 0.0, 259.0, 265.5, 265.5, 232.9, 249.6 ])) # Liquid density: rho = AD + BD*T + CD*T^2 ad: np.ndarray = field(default_factory=lambda: np.array([ 1.0, 1.0, 1.0, 23.3, 33.9, 32.8, 49.9, 50.5 ])) bd: np.ndarray = field(default_factory=lambda: np.array([ 0.0, 0.0, 0.0, -0.0700, -0.0957, -0.0995, -0.0191, -0.0541 ])) cd: np.ndarray = field(default_factory=lambda: np.array([ 0.0, 0.0, 0.0, -0.0002, -0.000152, -0.000233, -0.000425, -0.000150 ])) # Liquid enthalpy: H = AH*T + BH*T^2/2 + CH*T^3/3 ah: np.ndarray = field(default_factory=lambda: np.array([ 1.0e-6, 1.0e-6, 1.0e-6, 0.960e-6, 0.573e-6, 0.652e-6, 0.515e-6, 0.471e-6 ])) bh: np.ndarray = field(default_factory=lambda: np.array([ 0.0, 0.0, 0.0, 8.70e-9, 2.41e-9, 2.18e-9, 5.65e-10, 8.70e-10 ])) ch: np.ndarray = field(default_factory=lambda: np.array([ 0.0, 0.0, 0.0, 4.81e-11, 1.82e-11, 1.94e-11, 3.82e-12, 2.62e-12 ])) # Heat of vaporization av: np.ndarray = field(default_factory=lambda: np.array([ 1.0e-6, 1.0e-6, 1.0e-6, 86.7e-6, 160.0e-6, 160.0e-6, 225.0e-6, 209.0e-6 ])) # Gas heat capacity: Cp = AG + BG*T + CG*T^2 ag: np.ndarray = field(default_factory=lambda: np.array([ 3.411e-6, 0.3799e-6, 0.2491e-6, 0.3567e-6, 0.3463e-6, 0.3930e-6, 0.170e-6, 0.150e-6 ])) bg: np.ndarray = field(default_factory=lambda: np.array([ 7.18e-10, 1.08e-9, 1.36e-11, 8.51e-10, 8.96e-10, 1.02e-9, 0.0, 0.0 ])) cg: np.ndarray = field(default_factory=lambda: np.array([ 6.0e-13, -3.98e-13, -3.93e-14, -3.12e-13, -3.27e-13, -3.12e-13, 0.0, 0.0 ])) @dataclass class TEProcBlock: """Process state common block (/TEPROC/). Contains all internal process variables for reactor, separator, stripper (condenser), and compressor. """ # Reactor state uclr: np.ndarray = field(default_factory=lambda: np.zeros(8)) # Liquid component holdup ucvr: np.ndarray = field(default_factory=lambda: np.zeros(8)) # Vapor component holdup utlr: float = 0.0 # Total liquid holdup utvr: float = 0.0 # Total vapor holdup xlr: np.ndarray = field(default_factory=lambda: np.zeros(8)) # Liquid mole fractions xvr: np.ndarray = field(default_factory=lambda: np.zeros(8)) # Vapor mole fractions etr: float = 0.0 # Total energy esr: float = 0.0 # Specific energy tcr: float = 0.0 # Temperature (C) tkr: float = 0.0 # Temperature (K) dlr: float = 0.0 # Liquid density vlr: float = 0.0 # Liquid volume vvr: float = 0.0 # Vapor volume vtr: float = 1300.0 # Total volume ptr: float = 0.0 # Total pressure ppr: np.ndarray = field(default_factory=lambda: np.zeros(8)) # Partial pressures crxr: np.ndarray = field(default_factory=lambda: np.zeros(8)) # Reaction rates by component rr: np.ndarray = field(default_factory=lambda: np.zeros(4)) # Reaction rates rh: float = 0.0 # Heat of reaction fwr: float = 0.0 # Cooling water flow twr: float = 0.0 # Cooling water temperature qur: float = 0.0 # Heat transfer hwr: float = 7060.0 # Heat transfer coefficient uar: float = 0.0 # Overall heat transfer # Separator state ucls: np.ndarray = field(default_factory=lambda: np.zeros(8)) ucvs: np.ndarray = field(default_factory=lambda: np.zeros(8)) utls: float = 0.0 utvs: float = 0.0 xls: np.ndarray = field(default_factory=lambda: np.zeros(8)) xvs: np.ndarray = field(default_factory=lambda: np.zeros(8)) ets: float = 0.0 ess: float = 0.0 tcs: float = 0.0 tks: float = 0.0 dls: float = 0.0 vls: float = 0.0 vvs: float = 0.0 vts: float = 3500.0 pts: float = 0.0 pps: np.ndarray = field(default_factory=lambda: np.zeros(8)) fws: float = 0.0 tws: float = 0.0 qus: float = 0.0 hws: float = 11138.0 # Stripper/Condenser state uclc: np.ndarray = field(default_factory=lambda: np.zeros(8)) utlc: float = 0.0 xlc: np.ndarray = field(default_factory=lambda: np.zeros(8)) etc: float = 0.0 esc: float = 0.0 tcc: float = 0.0 dlc: float = 0.0 vlc: float = 0.0 vtc: float = 156.5 quc: float = 0.0 # Compressor state ucvv: np.ndarray = field(default_factory=lambda: np.zeros(8)) utvv: float = 0.0 xvv: np.ndarray = field(default_factory=lambda: np.zeros(8)) etv: float = 0.0 esv: float = 0.0 tcv: float = 0.0 tkv: float = 0.0 vtv: float = 5000.0 ptv: float = 0.0 # Valve states vcv: np.ndarray = field(default_factory=lambda: np.zeros(12)) # Valve command values vrng: np.ndarray = field(default_factory=lambda: np.array([ 400.0, 400.0, 100.0, 1500.0, 0.0, 0.0, 1500.0, 1000.0, 0.03, 1000.0, 1200.0, 0.0 ])) # Valve ranges vtau: np.ndarray = field(default_factory=lambda: np.array([ 8.0, 8.0, 6.0, 9.0, 7.0, 5.0, 5.0, 5.0, 120.0, 5.0, 5.0, 5.0 ]) / 3600.0) # Time constants (hours) # Stream flows and compositions ftm: np.ndarray = field(default_factory=lambda: np.zeros(13)) # Total molar flows fcm: np.ndarray = field(default_factory=lambda: np.zeros((8, 13))) # Component molar flows xst: np.ndarray = field(default_factory=lambda: np.zeros((8, 13))) # Stream compositions xmws: np.ndarray = field(default_factory=lambda: np.zeros(13)) # Stream molecular weights hst: np.ndarray = field(default_factory=lambda: np.zeros(13)) # Stream enthalpies tst: np.ndarray = field(default_factory=lambda: np.zeros(13)) # Stream temperatures sfr: np.ndarray = field(default_factory=lambda: np.zeros(8)) # Separation factors # Compressor parameters cpflmx: float = 280275.0 cpprmx: float = 1.3 cpdh: float = 0.0 # Cooling water inlet temperatures tcwr: float = 0.0 tcws: float = 0.0 # Heat of reaction and agitator htr: np.ndarray = field(default_factory=lambda: np.array([0.06899381054, 0.05, 0.0])) agsp: float = 0.0 # Measurement delays and noise xdel: np.ndarray = field(default_factory=lambda: np.zeros(41)) # Delayed measurements xns: np.ndarray = field(default_factory=lambda: np.zeros(41)) # Noise std devs tgas: float = 0.1 # Gas sampling time tprod: float = 0.25 # Product sampling time # Valve sticking vst: np.ndarray = field(default_factory=lambda: np.full(12, 2.0)) # Sticking threshold ivst: np.ndarray = field(default_factory=lambda: np.zeros(12, dtype=np.int32)) # Sticking flags @dataclass class WalkBlock: """Disturbance random walk common block (/WLK/). Contains parameters for cubic spline interpolation of random disturbances. """ adist: np.ndarray = field(default_factory=lambda: np.zeros(12)) bdist: np.ndarray = field(default_factory=lambda: np.zeros(12)) cdist: np.ndarray = field(default_factory=lambda: np.zeros(12)) ddist: np.ndarray = field(default_factory=lambda: np.zeros(12)) tlast: np.ndarray = field(default_factory=lambda: np.zeros(12)) tnext: np.ndarray = field(default_factory=lambda: np.full(12, 0.1)) hspan: np.ndarray = field(default_factory=lambda: np.array([ 0.2, 0.7, 0.25, 0.7, 0.15, 0.15, 1.0, 1.0, 0.4, 1.5, 2.0, 1.5 ])) hzero: np.ndarray = field(default_factory=lambda: np.array([ 0.5, 1.0, 0.5, 1.0, 0.25, 0.25, 2.0, 2.0, 0.5, 2.0, 3.0, 2.0 ])) sspan: np.ndarray = field(default_factory=lambda: np.array([ 0.03, 0.003, 10.0, 10.0, 10.0, 10.0, 0.25, 0.25, 0.25, 0.0, 0.0, 0.0 ])) szero: np.ndarray = field(default_factory=lambda: np.array([ 0.485, 0.005, 45.0, 45.0, 35.0, 40.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 ])) spspan: np.ndarray = field(default_factory=lambda: np.zeros(12)) idvwlk: np.ndarray = field(default_factory=lambda: np.zeros(12, dtype=np.int32)) # ============================================================================= # Pure Python TEP Process Implementation # ============================================================================= class PythonTEProcess: """Pure Python implementation of the Tennessee Eastman Process. This class replicates the original Fortran implementation (teprob.f) without any Fortran dependency. It provides the same interface as FortranTEProcess for drop-in replacement. The implementation includes: - Linear congruential random number generator (matching Fortran) - Full process dynamics for reactor, separator, stripper, compressor - Reaction kinetics and heat transfer - Valve dynamics with sticking - Measurement noise and sampling delays - 20 process disturbances Example: >>> process = PythonTEProcess(random_seed=1234) >>> process.initialize() >>> for _ in range(3600): # 1 hour ... process.step() >>> print(process.xmeas[:5]) """ def __init__(self, random_seed: Optional[int] = None): """Initialize the Python TEP process. Parameters ---------- random_seed : int, optional Random seed for reproducibility. If None, uses the default Fortran seed (4651207995). """ self._nn = 50 # Number of state variables self._initialized = False self.time = 0.0 self._random_seed = random_seed # Random number generator state self._g = 4651207995.0 # Default Fortran seed # Shutdown flag - once set, no further integration self._shutdown = False # State vectors self.yy = np.zeros(self._nn, dtype=np.float64) self.yp = np.zeros(self._nn, dtype=np.float64) # Common blocks self._const = ConstBlock() self._teproc = TEProcBlock() self._wlk = WalkBlock() # Process variables (PV common block) self._xmeas = np.zeros(41, dtype=np.float64) self._xmv = np.zeros(12, dtype=np.float64) # Disturbance vector (DVEC common block) self._idv = np.zeros(20, dtype=np.int32) # Create state wrapper for API compatibility self.state = PythonTEProcessState(self) # Create disturbances wrapper for API compatibility self.disturbances = PythonDisturbanceManager(self) # ========================================================================= # Random Number Generator (TESUB7) # ========================================================================= def _tesub7(self, i: int) -> float: """Linear congruential random number generator. Replicates Fortran TESUB7: G = MOD(G * 9228907, 2^32) if i >= 0: return G / 2^32 (range 0 to 1) if i < 0: return 2*G/2^32 - 1 (range -1 to 1) Parameters ---------- i : int If >= 0, returns value in [0, 1) If < 0, returns value in [-1, 1) Returns ------- float Random number """ # Use integer arithmetic to avoid float64 precision loss # The product can exceed 2^53, so we use Python's arbitrary precision int g_int = int(self._g) g_int = (g_int * 9228907) % 4294967296 self._g = float(g_int) if i >= 0: return self._g / 4294967296.0 else: return 2.0 * self._g / 4294967296.0 - 1.0 # ========================================================================= # Thermodynamic Functions (TESUB1-4) # ========================================================================= def _tesub1(self, z: np.ndarray, t: float, ity: int) -> float: """Calculate enthalpy from composition and temperature. Parameters ---------- z : np.ndarray Mole fractions (8 elements) t : float Temperature (deg C) ity : int 0 = liquid, 1 = gas, 2 = gas with pressure correction Returns ------- float Specific enthalpy """ const = self._const h = 0.0 if ity == 0: # Liquid enthalpy for i in range(8): hi = t * (const.ah[i] + const.bh[i] * t / 2.0 + const.ch[i] * t**2 / 3.0) hi = 1.8 * hi h += z[i] * const.xmw[i] * hi else: # Gas enthalpy for i in range(8): hi = t * (const.ag[i] + const.bg[i] * t / 2.0 + const.cg[i] * t**2 / 3.0) hi = 1.8 * hi hi += const.av[i] h += z[i] * const.xmw[i] * hi if ity == 2: # Pressure correction for compressor r = 3.57696e-6 h -= r * (t + 273.15) return h def _tesub2(self, z: np.ndarray, t_init: float, h: float, ity: int) -> float: """Calculate temperature from composition and enthalpy (Newton iteration). Parameters ---------- z : np.ndarray Mole fractions (8 elements) t_init : float Initial temperature guess (deg C) h : float Target specific enthalpy ity : int 0 = liquid, 1 = gas, 2 = gas with pressure correction Returns ------- float Temperature (deg C) """ t = t_init for _ in range(100): htest = self._tesub1(z, t, ity) err = htest - h dh = self._tesub3(z, t, ity) dt = -err / dh t += dt if abs(dt) < 1.0e-12: return t return t_init # Return initial guess if not converged def _tesub3(self, z: np.ndarray, t: float, ity: int) -> float: """Calculate enthalpy derivative dH/dT. Parameters ---------- z : np.ndarray Mole fractions (8 elements) t : float Temperature (deg C) ity : int 0 = liquid, 1 = gas, 2 = gas with pressure correction Returns ------- float dH/dT """ const = self._const dh = 0.0 if ity == 0: # Liquid for i in range(8): dhi = const.ah[i] + const.bh[i] * t + const.ch[i] * t**2 dhi = 1.8 * dhi dh += z[i] * const.xmw[i] * dhi else: # Gas for i in range(8): dhi = const.ag[i] + const.bg[i] * t + const.cg[i] * t**2 dhi = 1.8 * dhi dh += z[i] * const.xmw[i] * dhi if ity == 2: r = 3.57696e-6 dh -= r return dh def _tesub4(self, x: np.ndarray, t: float) -> float: """Calculate liquid density from composition and temperature. Parameters ---------- x : np.ndarray Mole fractions (8 elements) t : float Temperature (deg C) Returns ------- float Liquid density (kmol/m^3) """ const = self._const v = 0.0 for i in range(8): v += x[i] * const.xmw[i] / (const.ad[i] + (const.bd[i] + const.cd[i] * t) * t) return 1.0 / v if v > 0 else 1.0 # ========================================================================= # Disturbance Walk Functions (TESUB5, TESUB6, TESUB8) # ========================================================================= def _tesub5(self, s: float, sp: float, idx: int, idvflag: int): """Generate next random walk interval (cubic spline parameters). Updates the walk block parameters for disturbance index idx. Parameters ---------- s : float Current signal value sp : float Current signal derivative idx : int Disturbance index (0-based) idvflag : int Disturbance active flag (0 or 1) """ wlk = self._wlk # Generate random interval h = wlk.hspan[idx] * self._tesub7(-1) + wlk.hzero[idx] s1 = wlk.sspan[idx] * self._tesub7(-1) * idvflag + wlk.szero[idx] s1p = wlk.spspan[idx] * self._tesub7(-1) * idvflag # Calculate cubic spline coefficients wlk.adist[idx] = s wlk.bdist[idx] = sp wlk.cdist[idx] = (3.0 * (s1 - s) - h * (s1p + 2.0 * sp)) / h**2 wlk.ddist[idx] = (2.0 * (s - s1) + h * (s1p + sp)) / h**3 wlk.tnext[idx] = wlk.tlast[idx] + h def _tesub6(self, std: float) -> float: """Generate random noise with given standard deviation. Uses sum of 12 uniform random numbers (approximates normal distribution). Parameters ---------- std : float Standard deviation Returns ------- float Random noise value """ x = 0.0 for _ in range(12): x += self._tesub7(1) return (x - 6.0) * std def _tesub8(self, i: int, t: float) -> float: """Evaluate cubic spline for disturbance walk. Parameters ---------- i : int Disturbance index (1-based, as in Fortran) t : float Current time (hours) Returns ------- float Disturbance value """ wlk = self._wlk idx = i - 1 # Convert to 0-based h = t - wlk.tlast[idx] return wlk.adist[idx] + h * (wlk.bdist[idx] + h * (wlk.cdist[idx] + h * wlk.ddist[idx])) # ========================================================================= # Initialization (TEINIT) # ========================================================================= def _initialize(self): """Internal initialization method for TEPSimulator compatibility.""" self.initialize() def initialize(self): """Initialize the process to steady-state conditions. Replicates Fortran TEINIT subroutine. Sets initial state values, physical constants, and calls TEFUNC to compute initial derivatives. """ # Reset shutdown flag self._shutdown = False # Reset random seed self._g = 4651207995.0 # Initialize physical constants (already done in ConstBlock) # The constants are initialized in the dataclass defaults # Set initial state values (YY) self.yy = np.array([ 10.40491389, # YY(1) 4.363996017, # YY(2) 7.570059737, # YY(3) 0.4230042431, # YY(4) 24.15513437, # YY(5) 2.942597645, # YY(6) 154.3770655, # YY(7) 159.1865960, # YY(8) 2.808522723, # YY(9) 63.75581199, # YY(10) 26.74026066, # YY(11) 46.38532432, # YY(12) 0.2464521543, # YY(13) 15.20484404, # YY(14) 1.852266172, # YY(15) 52.44639459, # YY(16) 41.20394008, # YY(17) 0.5699317760, # YY(18) 0.4306056376, # YY(19) 7.9906200783e-03, # YY(20) 0.9056036089, # YY(21) 1.6054258216e-02, # YY(22) 0.7509759687, # YY(23) 8.8582855955e-02, # YY(24) 48.27726193, # YY(25) 39.38459028, # YY(26) 0.3755297257, # YY(27) 107.7562698, # YY(28) 29.77250546, # YY(29) 88.32481135, # YY(30) 23.03929507, # YY(31) 62.85848794, # YY(32) 5.546318688, # YY(33) 11.92244772, # YY(34) 5.555448243, # YY(35) 0.9218489762, # YY(36) 94.59927549, # YY(37) 77.29698353, # YY(38) 63.05263039, # YY(39) - XMV(1) 53.97970677, # YY(40) - XMV(2) 24.64355755, # YY(41) - XMV(3) 61.30192144, # YY(42) - XMV(4) 22.21000000, # YY(43) - XMV(5) 40.06374673, # YY(44) - XMV(6) 38.10034370, # YY(45) - XMV(7) 46.53415582, # YY(46) - XMV(8) 47.44573456, # YY(47) - XMV(9) 41.10581288, # YY(48) - XMV(10) 18.11349055, # YY(49) - XMV(11) 50.00000000, # YY(50) - XMV(12) ], dtype=np.float64) # Initialize XMV from state and valve parameters tp = self._teproc for i in range(12): self._xmv[i] = self.yy[38 + i] tp.vcv[i] = self._xmv[i] tp.vst[i] = 2.0 tp.ivst[i] = 0 # Initialize separation factors tp.sfr = np.array([0.99500, 0.99100, 0.99000, 0.91600, 0.93600, 0.93800, 0.05800, 0.03010]) # Initialize feed stream compositions # Stream 1: A Feed (mostly D) tp.xst[0, 0] = 0.0 tp.xst[1, 0] = 0.0001 tp.xst[2, 0] = 0.0 tp.xst[3, 0] = 0.9999 tp.xst[4, 0] = 0.0 tp.xst[5, 0] = 0.0 tp.xst[6, 0] = 0.0 tp.xst[7, 0] = 0.0 tp.tst[0] = 45.0 # Stream 2: D Feed (mostly E) tp.xst[0, 1] = 0.0 tp.xst[1, 1] = 0.0 tp.xst[2, 1] = 0.0 tp.xst[3, 1] = 0.0 tp.xst[4, 1] = 0.9999 tp.xst[5, 1] = 0.0001 tp.xst[6, 1] = 0.0 tp.xst[7, 1] = 0.0 tp.tst[1] = 45.0 # Stream 3: E Feed (mostly A) tp.xst[0, 2] = 0.9999 tp.xst[1, 2] = 0.0001 tp.xst[2, 2] = 0.0 tp.xst[3, 2] = 0.0 tp.xst[4, 2] = 0.0 tp.xst[5, 2] = 0.0 tp.xst[6, 2] = 0.0 tp.xst[7, 2] = 0.0 tp.tst[2] = 45.0 # Stream 4: A and C Feed tp.xst[0, 3] = 0.4850 tp.xst[1, 3] = 0.0050 tp.xst[2, 3] = 0.5100 tp.xst[3, 3] = 0.0 tp.xst[4, 3] = 0.0 tp.xst[5, 3] = 0.0 tp.xst[6, 3] = 0.0 tp.xst[7, 3] = 0.0 tp.tst[3] = 45.0 # Initialize measurement noise standard deviations tp.xns = np.array([ 0.0012, 18.000, 22.000, 0.0500, 0.2000, 0.2100, 0.3000, 0.5000, 0.0100, 0.0017, 0.0100, 1.0000, 0.3000, 0.1250, 1.0000, 0.3000, 0.1150, 0.0100, 1.1500, 0.2000, 0.0100, 0.0100, 0.250, 0.100, 0.250, 0.100, 0.250, 0.025, 0.250, 0.100, 0.250, 0.100, 0.250, 0.025, 0.050, 0.050, 0.010, 0.010, 0.010, 0.500, 0.500 ]) # Clear disturbances self._idv[:] = 0 # Initialize walk block wlk = self._wlk for i in range(12): wlk.tlast[i] = 0.0 wlk.tnext[i] = 0.1 wlk.adist[i] = wlk.szero[i] wlk.bdist[i] = 0.0 wlk.cdist[i] = 0.0 wlk.ddist[i] = 0.0 # Set time to zero self.time = 0.0 # Override random seed if provided if self._random_seed is not None: self._g = float(self._random_seed) # Call TEFUNC to compute initial derivatives self._tefunc() self._initialized = True # ========================================================================= # Main Simulation Function (TEFUNC) # ========================================================================= def _tefunc(self): """Compute state derivatives. Replicates Fortran TEFUNC subroutine. This is the core simulation function that computes all process dynamics. """ tp = self._teproc wlk = self._wlk const = self._const idv = self._idv time = self.time yy = self.yy yp = self.yp # Normalize IDV values to 0 or 1 for i in range(20): idv[i] = 1 if idv[i] > 0 else 0 # Map IDV to walk indices wlk.idvwlk[0] = idv[7] # IDV(8) wlk.idvwlk[1] = idv[7] # IDV(8) wlk.idvwlk[2] = idv[8] # IDV(9) wlk.idvwlk[3] = idv[9] # IDV(10) wlk.idvwlk[4] = idv[10] # IDV(11) wlk.idvwlk[5] = idv[11] # IDV(12) wlk.idvwlk[6] = idv[12] # IDV(13) wlk.idvwlk[7] = idv[12] # IDV(13) wlk.idvwlk[8] = idv[15] # IDV(16) wlk.idvwlk[9] = idv[16] # IDV(17) wlk.idvwlk[10] = idv[17] # IDV(18) wlk.idvwlk[11] = idv[19] # IDV(20) # Update random walks for disturbances 1-9 for i in range(9): if time >= wlk.tnext[i]: hwlk = wlk.tnext[i] - wlk.tlast[i] swlk = wlk.adist[i] + hwlk * (wlk.bdist[i] + hwlk * (wlk.cdist[i] + hwlk * wlk.ddist[i])) spwlk = wlk.bdist[i] + hwlk * (2.0 * wlk.cdist[i] + 3.0 * hwlk * wlk.ddist[i]) wlk.tlast[i] = wlk.tnext[i] self._tesub5(swlk, spwlk, i, wlk.idvwlk[i]) # Update random walks for disturbances 10-12 (special handling) for i in range(9, 12): if time >= wlk.tnext[i]: hwlk = wlk.tnext[i] - wlk.tlast[i] swlk = wlk.adist[i] + hwlk * (wlk.bdist[i] + hwlk * (wlk.cdist[i] + hwlk * wlk.ddist[i])) spwlk = wlk.bdist[i] + hwlk * (2.0 * wlk.cdist[i] + 3.0 * hwlk * wlk.ddist[i]) wlk.tlast[i] = wlk.tnext[i] if swlk > 0.1: wlk.adist[i] = swlk wlk.bdist[i] = spwlk wlk.cdist[i] = -(3.0 * swlk + 0.2 * spwlk) / 0.01 wlk.ddist[i] = (2.0 * swlk + 0.1 * spwlk) / 0.001 wlk.tnext[i] = wlk.tlast[i] + 0.1 else: hwlk = wlk.hspan[i] * self._tesub7(-1) + wlk.hzero[i] wlk.adist[i] = 0.0 wlk.bdist[i] = 0.0 wlk.cdist[i] = float(wlk.idvwlk[i]) / hwlk**2 wlk.ddist[i] = 0.0 wlk.tnext[i] = wlk.tlast[i] + hwlk # Reset walks at time = 0 if time == 0.0: for i in range(12): wlk.adist[i] = wlk.szero[i] wlk.bdist[i] = 0.0 wlk.cdist[i] = 0.0 wlk.ddist[i] = 0.0 wlk.tlast[i] = 0.0 wlk.tnext[i] = 0.1 # Apply disturbances to feed compositions tp.xst[0, 3] = self._tesub8(1, time) - idv[0] * 0.03 - idv[1] * 2.43719e-3 tp.xst[1, 3] = self._tesub8(2, time) + idv[1] * 0.005 tp.xst[2, 3] = 1.0 - tp.xst[0, 3] - tp.xst[1, 3] tp.tst[0] = self._tesub8(3, time) + idv[2] * 5.0 tp.tst[3] = self._tesub8(4, time) tp.tcwr = self._tesub8(5, time) + idv[3] * 5.0 tp.tcws = self._tesub8(6, time) + idv[4] * 5.0 r1f = self._tesub8(7, time) r2f = self._tesub8(8, time) # Extract state variables # Reactor vapor holdup (components 1-3) for i in range(3): tp.ucvr[i] = yy[i] tp.ucvs[i] = yy[i + 9] tp.uclr[i] = 0.0 tp.ucls[i] = 0.0 # Reactor liquid holdup (components 4-8) for i in range(3, 8): tp.uclr[i] = yy[i] tp.ucls[i] = yy[i + 9] # Condenser and compressor holdup for i in range(8): tp.uclc[i] = yy[i + 18] tp.ucvv[i] = yy[i + 27] # Energy states tp.etr = yy[8] tp.ets = yy[17] tp.etc = yy[26] tp.etv = yy[35] # Cooling water temperatures tp.twr = yy[36] tp.tws = yy[37] # Valve positions vpos = yy[38:50].copy() # Calculate total holdups tp.utlr = np.sum(tp.uclr) tp.utls = np.sum(tp.ucls) tp.utlc = np.sum(tp.uclc) tp.utvv = np.sum(tp.ucvv) # Calculate mole fractions for i in range(8): tp.xlr[i] = tp.uclr[i] / tp.utlr if tp.utlr > 0 else 0.0 tp.xls[i] = tp.ucls[i] / tp.utls if tp.utls > 0 else 0.0 tp.xlc[i] = tp.uclc[i] / tp.utlc if tp.utlc > 0 else 0.0 tp.xvv[i] = tp.ucvv[i] / tp.utvv if tp.utvv > 0 else 0.0 # Calculate specific energies tp.esr = tp.etr / tp.utlr if tp.utlr > 0 else 0.0 tp.ess = tp.ets / tp.utls if tp.utls > 0 else 0.0 tp.esc = tp.etc / tp.utlc if tp.utlc > 0 else 0.0 tp.esv = tp.etv / tp.utvv if tp.utvv > 0 else 0.0 # Calculate temperatures from enthalpies tp.tcr = self._tesub2(tp.xlr, tp.tcr if tp.tcr > 0 else 120.0, tp.esr, 0) tp.tkr = tp.tcr + 273.15 tp.tcs = self._tesub2(tp.xls, tp.tcs if tp.tcs > 0 else 80.0, tp.ess, 0) tp.tks = tp.tcs + 273.15 tp.tcc = self._tesub2(tp.xlc, tp.tcc if tp.tcc > 0 else 65.0, tp.esc, 0) tp.tcv = self._tesub2(tp.xvv, tp.tcv if tp.tcv > 0 else 100.0, tp.esv, 2) tp.tkv = tp.tcv + 273.15 # Calculate densities tp.dlr = self._tesub4(tp.xlr, tp.tcr) tp.dls = self._tesub4(tp.xls, tp.tcs) tp.dlc = self._tesub4(tp.xlc, tp.tcc) # Calculate volumes tp.vlr = tp.utlr / tp.dlr if tp.dlr > 0 else 0.0 tp.vls = tp.utls / tp.dls if tp.dls > 0 else 0.0 tp.vlc = tp.utlc / tp.dlc if tp.dlc > 0 else 0.0 tp.vvr = tp.vtr - tp.vlr tp.vvs = tp.vts - tp.vls # Gas constant rg = 998.9 # Calculate pressures (reactor and separator) tp.ptr = 0.0 tp.pts = 0.0 # Non-condensable components (ideal gas) for i in range(3): tp.ppr[i] = tp.ucvr[i] * rg * tp.tkr / tp.vvr if tp.vvr > 0 else 0.0 tp.ptr += tp.ppr[i] tp.pps[i] = tp.ucvs[i] * rg * tp.tks / tp.vvs if tp.vvs > 0 else 0.0 tp.pts += tp.pps[i] # Condensable components (vapor pressure) for i in range(3, 8): vpr = np.exp(const.avp[i] + const.bvp[i] / (tp.tcr + const.cvp[i])) tp.ppr[i] = vpr * tp.xlr[i] tp.ptr += tp.ppr[i] vpr = np.exp(const.avp[i] + const.bvp[i] / (tp.tcs + const.cvp[i])) tp.pps[i] = vpr * tp.xls[i] tp.pts += tp.pps[i] # Compressor pressure tp.ptv = tp.utvv * rg * tp.tkv / tp.vtv if tp.vtv > 0 else 0.0 # Calculate vapor compositions for i in range(8): tp.xvr[i] = tp.ppr[i] / tp.ptr if tp.ptr > 0 else 0.0 tp.xvs[i] = tp.pps[i] / tp.pts if tp.pts > 0 else 0.0 # Calculate total vapor holdups tp.utvr = tp.ptr * tp.vvr / rg / tp.tkr if (rg * tp.tkr) > 0 else 0.0 tp.utvs = tp.pts * tp.vvs / rg / tp.tks if (rg * tp.tks) > 0 else 0.0 # Update condensable vapor holdups for i in range(3, 8): tp.ucvr[i] = tp.utvr * tp.xvr[i] tp.ucvs[i] = tp.utvs * tp.xvs[i] # Reaction kinetics tp.rr[0] = np.exp(31.5859536 - 40000.0 / 1.987 / tp.tkr) * r1f tp.rr[1] = np.exp(3.00094014 - 20000.0 / 1.987 / tp.tkr) * r2f tp.rr[2] = np.exp(53.4060443 - 60000.0 / 1.987 / tp.tkr) tp.rr[3] = tp.rr[2] * 0.767488334 if tp.ppr[0] > 0.0 and tp.ppr[2] > 0.0: r1f_pp = tp.ppr[0] ** 1.1544 r2f_pp = tp.ppr[2] ** 0.3735 tp.rr[0] = tp.rr[0] * r1f_pp * r2f_pp * tp.ppr[3] tp.rr[1] = tp.rr[1] * r1f_pp * r2f_pp * tp.ppr[4] else: tp.rr[0] = 0.0 tp.rr[1] = 0.0 tp.rr[2] = tp.rr[2] * tp.ppr[0] * tp.ppr[4] tp.rr[3] = tp.rr[3] * tp.ppr[0] * tp.ppr[3] # Scale by vapor volume for i in range(4): tp.rr[i] = tp.rr[i] * tp.vvr # Component reaction rates tp.crxr[0] = -tp.rr[0] - tp.rr[1] - tp.rr[2] tp.crxr[1] = 0.0 tp.crxr[2] = -tp.rr[0] - tp.rr[1] tp.crxr[3] = -tp.rr[0] - 1.5 * tp.rr[3] tp.crxr[4] = -tp.rr[1] - tp.rr[2] tp.crxr[5] = tp.rr[2] + tp.rr[3] tp.crxr[6] = tp.rr[0] tp.crxr[7] = tp.rr[1] # Heat of reaction tp.rh = tp.rr[0] * tp.htr[0] + tp.rr[1] * tp.htr[1] # Stream compositions and molecular weights tp.xmws[0] = 0.0 tp.xmws[1] = 0.0 tp.xmws[5] = 0.0 tp.xmws[7] = 0.0 tp.xmws[8] = 0.0 tp.xmws[9] = 0.0 for i in range(8): tp.xst[i, 5] = tp.xvv[i] tp.xst[i, 7] = tp.xvr[i] tp.xst[i, 8] = tp.xvs[i] tp.xst[i, 9] = tp.xvs[i] tp.xst[i, 10] = tp.xls[i] tp.xst[i, 12] = tp.xlc[i] tp.xmws[0] += tp.xst[i, 0] * const.xmw[i] tp.xmws[1] += tp.xst[i, 1] * const.xmw[i] tp.xmws[5] += tp.xst[i, 5] * const.xmw[i] tp.xmws[7] += tp.xst[i, 7] * const.xmw[i] tp.xmws[8] += tp.xst[i, 8] * const.xmw[i] tp.xmws[9] += tp.xst[i, 9] * const.xmw[i] # Stream temperatures tp.tst[5] = tp.tcv tp.tst[7] = tp.tcr tp.tst[8] = tp.tcs tp.tst[9] = tp.tcs tp.tst[10] = tp.tcs tp.tst[12] = tp.tcc # Calculate stream enthalpies tp.hst[0] = self._tesub1(tp.xst[:, 0], tp.tst[0], 1) tp.hst[1] = self._tesub1(tp.xst[:, 1], tp.tst[1], 1) tp.hst[2] = self._tesub1(tp.xst[:, 2], tp.tst[2], 1) tp.hst[3] = self._tesub1(tp.xst[:, 3], tp.tst[3], 1) tp.hst[5] = self._tesub1(tp.xst[:, 5], tp.tst[5], 1) tp.hst[7] = self._tesub1(tp.xst[:, 7], tp.tst[7], 1) tp.hst[8] = self._tesub1(tp.xst[:, 8], tp.tst[8], 1) tp.hst[9] = tp.hst[8] tp.hst[10] = self._tesub1(tp.xst[:, 10], tp.tst[10], 0) tp.hst[12] = self._tesub1(tp.xst[:, 12], tp.tst[12], 0) # Calculate flows tp.ftm[0] = vpos[0] * tp.vrng[0] / 100.0 tp.ftm[1] = vpos[1] * tp.vrng[1] / 100.0 tp.ftm[2] = vpos[2] * (1.0 - idv[5]) * tp.vrng[2] / 100.0 tp.ftm[3] = vpos[3] * (1.0 - idv[6] * 0.2) * tp.vrng[3] / 100.0 + 1.0e-10 tp.ftm[10] = vpos[6] * tp.vrng[6] / 100.0 tp.ftm[12] = vpos[7] * tp.vrng[7] / 100.0 uac = vpos[8] * tp.vrng[8] * (1.0 + self._tesub8(9, time)) / 100.0 tp.fwr = vpos[9] * tp.vrng[9] / 100.0 tp.fws = vpos[10] * tp.vrng[10] / 100.0 tp.agsp = (vpos[11] + 150.0) / 100.0 # Pressure-driven flows dlp = tp.ptv - tp.ptr if dlp < 0.0: dlp = 0.0 flms = 1937.6 * np.sqrt(dlp) tp.ftm[5] = flms / tp.xmws[5] if tp.xmws[5] > 0 else 0.0 dlp = tp.ptr - tp.pts if dlp < 0.0: dlp = 0.0 flms = 4574.21 * np.sqrt(dlp) * (1.0 - 0.25 * self._tesub8(12, time)) tp.ftm[7] = flms / tp.xmws[7] if tp.xmws[7] > 0 else 0.0 dlp = tp.pts - 760.0 if dlp < 0.0: dlp = 0.0 flms = vpos[5] * 0.151169 * np.sqrt(dlp) tp.ftm[9] = flms / tp.xmws[9] if tp.xmws[9] > 0 else 0.0 # Compressor pr = tp.ptv / tp.pts if tp.pts > 0 else 1.0 if pr < 1.0: pr = 1.0 if pr > tp.cpprmx: pr = tp.cpprmx flcoef = tp.cpflmx / 1.197 flms = tp.cpflmx + flcoef * (1.0 - pr**3) tp.cpdh = flms * (tp.tcs + 273.15) * 1.8e-6 * 1.9872 * (tp.ptv - tp.pts) / (tp.xmws[8] * tp.pts) if (tp.xmws[8] * tp.pts) > 0 else 0.0 dlp = tp.ptv - tp.pts if dlp < 0.0: dlp = 0.0 flms = flms - vpos[4] * 53.349 * np.sqrt(dlp) if flms < 1.0e-3: flms = 1.0e-3 tp.ftm[8] = flms / tp.xmws[8] if tp.xmws[8] > 0 else 0.0 tp.hst[8] = tp.hst[8] + tp.cpdh / tp.ftm[8] if tp.ftm[8] > 0 else tp.hst[8] # Component flows for i in range(8): tp.fcm[i, 0] = tp.xst[i, 0] * tp.ftm[0] tp.fcm[i, 1] = tp.xst[i, 1] * tp.ftm[1] tp.fcm[i, 2] = tp.xst[i, 2] * tp.ftm[2] tp.fcm[i, 3] = tp.xst[i, 3] * tp.ftm[3] tp.fcm[i, 5] = tp.xst[i, 5] * tp.ftm[5] tp.fcm[i, 7] = tp.xst[i, 7] * tp.ftm[7] tp.fcm[i, 8] = tp.xst[i, 8] * tp.ftm[8] tp.fcm[i, 9] = tp.xst[i, 9] * tp.ftm[9] tp.fcm[i, 10] = tp.xst[i, 10] * tp.ftm[10] tp.fcm[i, 12] = tp.xst[i, 12] * tp.ftm[12] # Stripper separation if tp.ftm[10] > 0.1: if tp.tcc > 170.0: tmpfac = tp.tcc - 120.262 elif tp.tcc < 5.292: tmpfac = 0.1 else: tmpfac = 363.744 / (177.0 - tp.tcc) - 2.22579488 vovrl = tp.ftm[3] / tp.ftm[10] * tmpfac tp.sfr[3] = 8.5010 * vovrl / (1.0 + 8.5010 * vovrl) tp.sfr[4] = 11.402 * vovrl / (1.0 + 11.402 * vovrl) tp.sfr[5] = 11.795 * vovrl / (1.0 + 11.795 * vovrl) tp.sfr[6] = 0.0480 * vovrl / (1.0 + 0.0480 * vovrl) tp.sfr[7] = 0.0242 * vovrl / (1.0 + 0.0242 * vovrl) else: tp.sfr[3] = 0.9999 tp.sfr[4] = 0.999 tp.sfr[5] = 0.999 tp.sfr[6] = 0.99 tp.sfr[7] = 0.98 # Stripper inlet flows fin = np.zeros(8) for i in range(8): fin[i] = tp.fcm[i, 3] + tp.fcm[i, 10] # Stripper separation tp.ftm[4] = 0.0 tp.ftm[11] = 0.0 for i in range(8): tp.fcm[i, 4] = tp.sfr[i] * fin[i] tp.fcm[i, 11] = fin[i] - tp.fcm[i, 4] tp.ftm[4] += tp.fcm[i, 4] tp.ftm[11] += tp.fcm[i, 11] # Stream compositions for i in range(8): tp.xst[i, 4] = tp.fcm[i, 4] / tp.ftm[4] if tp.ftm[4] > 0 else 0.0 tp.xst[i, 11] = tp.fcm[i, 11] / tp.ftm[11] if tp.ftm[11] > 0 else 0.0 tp.tst[4] = tp.tcc tp.tst[11] = tp.tcc tp.hst[4] = self._tesub1(tp.xst[:, 4], tp.tst[4], 1) tp.hst[11] = self._tesub1(tp.xst[:, 11], tp.tst[11], 0) # Stream 7 = Stream 6 tp.ftm[6] = tp.ftm[5] tp.hst[6] = tp.hst[5] tp.tst[6] = tp.tst[5] for i in range(8): tp.xst[i, 6] = tp.xst[i, 5] tp.fcm[i, 6] = tp.fcm[i, 5] # Heat transfer calculations if tp.vlr / 7.8 > 50.0: uarlev = 1.0 elif tp.vlr / 7.8 < 10.0: uarlev = 0.0 else: uarlev = 0.025 * tp.vlr / 7.8 - 0.25 tp.uar = uarlev * (-0.5 * tp.agsp**2 + 2.75 * tp.agsp - 2.5) * 855490.0e-6 tp.qur = tp.uar * (tp.twr - tp.tcr) * (1.0 - 0.35 * self._tesub8(10, time)) uas = 0.404655 * (1.0 - 1.0 / (1.0 + (tp.ftm[7] / 3528.73)**4)) tp.qus = uas * (tp.tws - tp.tst[7]) * (1.0 - 0.25 * self._tesub8(11, time)) tp.quc = 0.0 if tp.tcc < 100.0: tp.quc = uac * (100.0 - tp.tcc) # Calculate measurements (without noise first) self._xmeas[0] = tp.ftm[2] * 0.359 / 35.3145 self._xmeas[1] = tp.ftm[0] * tp.xmws[0] * 0.454 self._xmeas[2] = tp.ftm[1] * tp.xmws[1] * 0.454 self._xmeas[3] = tp.ftm[3] * 0.359 / 35.3145 self._xmeas[4] = tp.ftm[8] * 0.359 / 35.3145 self._xmeas[5] = tp.ftm[5] * 0.359 / 35.3145 self._xmeas[6] = (tp.ptr - 760.0) / 760.0 * 101.325 self._xmeas[7] = (tp.vlr - 84.6) / 666.7 * 100.0 self._xmeas[8] = tp.tcr self._xmeas[9] = tp.ftm[9] * 0.359 / 35.3145 self._xmeas[10] = tp.tcs self._xmeas[11] = (tp.vls - 27.5) / 290.0 * 100.0 self._xmeas[12] = (tp.pts - 760.0) / 760.0 * 101.325 self._xmeas[13] = tp.ftm[10] / tp.dls / 35.3145 if tp.dls > 0 else 0.0 self._xmeas[14] = (tp.vlc - 78.25) / tp.vtc * 100.0 self._xmeas[15] = (tp.ptv - 760.0) / 760.0 * 101.325 self._xmeas[16] = tp.ftm[12] / tp.dlc / 35.3145 if tp.dlc > 0 else 0.0 self._xmeas[17] = tp.tcc self._xmeas[18] = tp.quc * 1.04e3 * 0.454 self._xmeas[19] = tp.cpdh * 0.29307e3 self._xmeas[20] = tp.twr self._xmeas[21] = tp.tws # Safety shutdown check isd = 0 if self._xmeas[6] > 3000.0: isd = 1 if tp.vlr / 35.3145 > 24.0: isd = 1 if tp.vlr / 35.3145 < 2.0: isd = 1 if self._xmeas[8] > 175.0: isd = 1 if tp.vls / 35.3145 > 12.0: isd = 1 if tp.vls / 35.3145 < 1.0: isd = 1 if tp.vlc / 35.3145 > 8.0: isd = 1 if tp.vlc / 35.3145 < 1.0: isd = 1 # Add measurement noise (only after time 0) if time > 0.0 and isd == 0: for i in range(22): self._xmeas[i] += self._tesub6(tp.xns[i]) # Sampled composition measurements (with delay) xcmp = np.zeros(41) xcmp[22] = tp.xst[0, 6] * 100.0 xcmp[23] = tp.xst[1, 6] * 100.0 xcmp[24] = tp.xst[2, 6] * 100.0 xcmp[25] = tp.xst[3, 6] * 100.0 xcmp[26] = tp.xst[4, 6] * 100.0 xcmp[27] = tp.xst[5, 6] * 100.0 xcmp[28] = tp.xst[0, 9] * 100.0 xcmp[29] = tp.xst[1, 9] * 100.0 xcmp[30] = tp.xst[2, 9] * 100.0 xcmp[31] = tp.xst[3, 9] * 100.0 xcmp[32] = tp.xst[4, 9] * 100.0 xcmp[33] = tp.xst[5, 9] * 100.0 xcmp[34] = tp.xst[6, 9] * 100.0 xcmp[35] = tp.xst[7, 9] * 100.0 xcmp[36] = tp.xst[3, 12] * 100.0 xcmp[37] = tp.xst[4, 12] * 100.0 xcmp[38] = tp.xst[5, 12] * 100.0 xcmp[39] = tp.xst[6, 12] * 100.0 xcmp[40] = tp.xst[7, 12] * 100.0 if time == 0.0: for i in range(22, 41): tp.xdel[i] = xcmp[i] self._xmeas[i] = xcmp[i] tp.tgas = 0.1 tp.tprod = 0.25 if time >= tp.tgas: for i in range(22, 36): self._xmeas[i] = tp.xdel[i] + self._tesub6(tp.xns[i]) tp.xdel[i] = xcmp[i] tp.tgas += 0.1 if time >= tp.tprod: for i in range(36, 41): self._xmeas[i] = tp.xdel[i] + self._tesub6(tp.xns[i]) tp.xdel[i] = xcmp[i] tp.tprod += 0.25 # State derivatives # Reactor component balances for i in range(8): yp[i] = tp.fcm[i, 6] - tp.fcm[i, 7] + tp.crxr[i] yp[i + 9] = tp.fcm[i, 7] - tp.fcm[i, 8] - tp.fcm[i, 9] - tp.fcm[i, 10] yp[i + 18] = tp.fcm[i, 11] - tp.fcm[i, 12] yp[i + 27] = tp.fcm[i, 0] + tp.fcm[i, 1] + tp.fcm[i, 2] + tp.fcm[i, 4] + tp.fcm[i, 8] - tp.fcm[i, 5] # Energy balances yp[8] = tp.hst[6] * tp.ftm[6] - tp.hst[7] * tp.ftm[7] + tp.rh + tp.qur yp[17] = tp.hst[7] * tp.ftm[7] - tp.hst[8] * tp.ftm[8] - tp.hst[9] * tp.ftm[9] - tp.hst[10] * tp.ftm[10] + tp.qus yp[26] = tp.hst[3] * tp.ftm[3] + tp.hst[10] * tp.ftm[10] - tp.hst[4] * tp.ftm[4] - tp.hst[12] * tp.ftm[12] + tp.quc yp[35] = tp.hst[0] * tp.ftm[0] + tp.hst[1] * tp.ftm[1] + tp.hst[2] * tp.ftm[2] + tp.hst[4] * tp.ftm[4] + tp.hst[8] * tp.ftm[8] - tp.hst[5] * tp.ftm[5] # Cooling water temperatures yp[36] = (tp.fwr * 500.53 * (tp.tcwr - tp.twr) - tp.qur * 1.0e6 / 1.8) / tp.hwr yp[37] = (tp.fws * 500.53 * (tp.tcws - tp.tws) - tp.qus * 1.0e6 / 1.8) / tp.hws # Valve sticking tp.ivst[9] = idv[13] # IDV(14) tp.ivst[10] = idv[14] # IDV(15) tp.ivst[4] = idv[18] # IDV(19) tp.ivst[6] = idv[18] # IDV(19) tp.ivst[7] = idv[18] # IDV(19) tp.ivst[8] = idv[18] # IDV(19) # Valve dynamics for i in range(12): if time == 0.0 or abs(tp.vcv[i] - self._xmv[i]) > tp.vst[i] * tp.ivst[i]: tp.vcv[i] = self._xmv[i] tp.vcv[i] = np.clip(tp.vcv[i], 0.0, 100.0) yp[38 + i] = (tp.vcv[i] - vpos[i]) / tp.vtau[i] # Shutdown: zero all derivatives and set flag if isd != 0: yp[:] = 0.0 self._shutdown = True # ========================================================================= # Public Interface Methods # ========================================================================= def step(self, dt: float = 1.0/3600.0): """Single Euler integration step. Parameters ---------- dt : float Time step in hours. Default is 1 second (1/3600 hours). """ if not self._initialized: raise RuntimeError("Process not initialized. Call initialize() first.") # Don't integrate if already shutdown if self._shutdown: self.time += dt return # Call tefunc to get derivatives self._tefunc() # If shutdown was triggered in tefunc, don't integrate if self._shutdown: self.time += dt return # Check for numerical instability before integration if not np.all(np.isfinite(self.yp)): self._shutdown = True self.yp[:] = 0.0 self.time += dt return # Euler integration: yy = yy + yp * dt self.yy[:] = self.yy + self.yp * dt self.time += dt # Check for numerical instability after integration if not np.all(np.isfinite(self.yy)): self._shutdown = True return # Apply valve constraints self._apply_constraints() def _apply_constraints(self): """Apply valve constraints (0-100%).""" for i in range(11): self._xmv[i] = np.clip(self._xmv[i], 0.0, 100.0) @property def xmeas(self) -> np.ndarray: """Get current measurement values.""" return self._xmeas.copy() @property def xmv(self) -> np.ndarray: """Get current manipulated variables.""" return self._xmv.copy() @property def idv(self) -> np.ndarray: """Get current disturbance vector.""" return self._idv.copy() def set_xmv(self, index: int, value: float): """Set a manipulated variable. Parameters ---------- index : int 1-based index of the manipulated variable (1-12). value : float Value to set (will be clipped to 0-100). """ if not 1 <= index <= 12: raise ValueError(f"XMV index must be 1-12, got {index}") self._xmv[index - 1] = np.clip(value, 0.0, 100.0) def set_idv(self, index: int, value: int): """Set a disturbance variable. Parameters ---------- index : int 1-based index of the disturbance (1-20). value : int 0 to disable, 1 to enable the disturbance. """ if not 1 <= index <= 20: raise ValueError(f"IDV index must be 1-20, got {index}") self._idv[index - 1] = value def clear_disturbances(self): """Clear all disturbances (set IDV to 0).""" self._idv[:] = 0 def get_xmeas(self) -> np.ndarray: """Get current measurement values (for TEPSimulator compatibility).""" return self._xmeas.copy() def get_xmv(self) -> np.ndarray: """Get current manipulated variables (for TEPSimulator compatibility).""" return self._xmv.copy() def evaluate(self, time: float, yy: np.ndarray) -> np.ndarray: """Evaluate derivatives using TEFUNC. Parameters ---------- time : float Current time in hours. yy : np.ndarray Current state vector. Returns ------- np.ndarray Derivative vector. """ self.yy[:] = yy self.time = time self._tefunc() return self.yp.copy() def is_shutdown(self) -> bool: """Check if process is in shutdown state.""" # Return stored flag first (handles numerical instability) if self._shutdown: return True # Also check shutdown conditions from measurements tp = self._teproc if self._xmeas[6] > 3000.0: return True if tp.vlr / 35.3145 > 24.0 or tp.vlr / 35.3145 < 2.0: return True if self._xmeas[8] > 175.0: return True if tp.vls / 35.3145 > 12.0 or tp.vls / 35.3145 < 1.0: return True if tp.vlc / 35.3145 > 8.0 or tp.vlc / 35.3145 < 1.0: return True return False def get_state(self) -> np.ndarray: """Get current state vector.""" return self.yy.copy() def set_state(self, state: np.ndarray): """Set state vector.""" if len(state) != self._nn: raise ValueError(f"State must have {self._nn} elements") self.yy[:] = state class PythonTEProcessState: """State wrapper that mimics TEProcess state interface.""" def __init__(self, process: 'PythonTEProcess'): self._process = process @property def yy(self) -> np.ndarray: return self._process.yy @yy.setter def yy(self, value: np.ndarray): self._process.yy[:] = value @property def xmeas(self) -> np.ndarray: return self._process.xmeas @property def xmv(self) -> np.ndarray: return self._process.xmv class PythonDisturbanceManager: """Disturbance manager wrapper for API compatibility.""" def __init__(self, process: 'PythonTEProcess'): self._process = process def clear_all_disturbances(self): """Clear all disturbances (set IDV to 0).""" self._process.clear_disturbances() def set_disturbance(self, index: int, value: int): """Set a disturbance variable.""" self._process.set_idv(index, value)