PIDretuner for Integrating Processes

  ! 340 total lines, October 13, 2025. 

  I M P L (c)

  Copyright and Property of i n d u s t r i @ l g o r i t h m s

&sSetting,@rValue
WRITECALCS,1
WRITEDATA,0

NUMLAGS,50
&sSetting,@rValue

&sMemory,@rValue
LENSSVAL,50000000
&sMemory,@rValue

! The PID data-vector (pid[]) with lower and upper ranges and step or stride increments along with the upper norm type and limit
! (unt, unl) and the SISOARX denominator degree and numerator degree and delay lowers and uppers (ml, mu, nl, nu, kl and ku).

Include-@sFile_Name
INCLUDEFILE1,DATAFRAME
Include-@sFile_Name

!Process plus PID closed-loop feedback data i.e., sp[], pv[], op[] (spo tuple) data-vectors in any arrangement, permutation or order.

Include-@sFile_Name
INCLUDEFILE2,DATAFRAME
Include-@sFile_Name
 
&sCalc,@sValue

! RETUNEPID()'s upper norm type i.e., (1) 1-norm, (2, default) 2-norm or (3) oo-norm.
unt,pid[22] !1, 2 or 3
unt,0*LT(unt;2)+5*EQ(unt;2)+10*GT(unt;2) ! Converted to a index or indice.

! RETUNEPID()'s upper norm limit (unl) (formulas flag) for the RETUNEPID()'s input-move-error (IME) upper norm limit for 
! self-regulating processes or output-error (OE) and negative (-ve) for integrating processes.
unl,pid[23] !-1000.0

! SISOARX() structure tuples where "m" is the denominator degree (0..oo), "n" is the numerator degree (0..oo) and "k" is the delay or 
! dead-time which must have a minimum value of one (1..oo) considering the zero-order-hold and discrete sampling.
!
! * Note that the dead-time or delay-time index (k) can be set to zero (0), but this is unrealistic as it implies that the process
!   output (y,t) is immediately and directly influenced by the process input (u,t) which is physically impossible especially in
!   discrete-time sampled systems.
ml,pid[24] !0
mu,pid[25] !0
nl,pid[26] !0
nu,pid[27] !0
kl,pid[28] !1
ku,pid[29] !10

! Calculate the length, size or cardinality of the ARX structure ranges or domains.

nm,mu-ml+1
nn,nu-nl+1
nk,ku-kl+1
&sCalc,@sValue

! Copy the SPO (sp,pv,op) tuple to the RYU (r,y,u) (mean-, steady-state- or opening-value-centered) tuple of data-vectors.

&sDataData,@sValue
r,sp
y,pv
u,op
&sDataData,@sValue

! Compute the means and norms of output- and input-move-error (oe, ime) performance.

&sDataCalc,@sValue
nt,LENGTH(r)

! * Note that for integrating processes, we require the opening, initial, starting or beginning value instead of the 
!   mean or steady-state value.

meanr,r[1]
meany,y[1]
meanu,u[1]

e1,NORMOE(r;y;1)
e2,NORMOE(r;y;2)
eoo,NORMOE(r;y;3)
du1,NORMIME(u;1)
du2,NORMIME(u;2)
duoo,NORMIME(u;3)
&sDataCalc,@sValue

&sDataData,@sValue

! Mean-center the process data into perturbation / deviation form by subtracting its mean, steady-state or 
! opening / initial / starting / beginning value.
 
r,r-meanr
y,y-meany
u,u-meanu

! Adjust the lower and upper process limits in the PID data-vector for the input and output by its mean, 
! steady-state or opening / initial / starting / beginning value.

pid,SWAP(pid[18]-meanu;18)
pid,SWAP(pid[19]-meanu;19)
pid,SWAP(pid[20]-meany;20)
pid,SWAP(pid[21]-meany;21)
&sDataData,@sValue

! Following the integrating process modeling of Kelly, J.D., "Tuning Digital PI Controllers for Minimal Variance in 
! Manipulated Input Moves applied to Imbalanced Systems with Delay", Canadian Journal of Chemical Engineering, 76, 
! (1998), we difference the integrating process output by two (2) or twice and by one (1) for the process input. 
!
! A discussion on the reasoning behind double differencing the output (level, height, holdup, accumulation, etc.) is
! given in Kelly (1998) as is based on assuming a "random walk" or IMA(1,0) for either the flow in or flow out depending
! on which is the unmeasured disturbance variable. 

&sDataData,@sValue
yd2,DIFFERENCE(y;2)
ud1,DIFFERENCE(u;1)
&sDataData,@sValue

! Automate the structure selection with SISOARX() using a simple brute-force enumerative search.
!
! SISO ARX lagged parametric model v. the transfer function (rational) model i.e.,
! y,t = y,t+0 = y,t-0 = a1 * y,t-1 + a2 * y,t-2 + ... + am * y,t-m + b0 * u,t-0-k + ... + bn+1 * u,t-n-k + e,t
!
! However, for integrating or imbalanced systems, Kelly (1998) formulates the level, holdup or pressure process system 
! as defined below:
!
! h,t = 2 * h,t-1 - 1 * h,t-2 + w0 * qo,t-k - w0 * qo,t-k-1 + a,t 
!
! or
!
! D^2 * h,t = w0 * D * qo,t-k + a,t 
!
! or
!
! D * h,t = w0 * qo,t-k + D^-1 * a,t
!
! where h,t is the height, level or pressure, qo,t is the outlet flow (or qi,t for inlet flow), D is the difference operator 
! i.e., D = 1-z^-1 and w0 is the process gain "omega" in the Box-Jenkins (BXJK) nomenclature.
!
! From Kelly (1998), the a,t time-series equals -w0 * D * qi,t-1 therefore the flow into the vessel unit is qi,t = qi,t-1 + a,t / -w0 
! where -w0 = tpd / area when the level or height is in length unit-of-measures.
!
! For interest, Kelly (1998) assumed that the inlet flow disturbance is represented as a "random walk" time-series which is denoted 
! as qi,t-1 = a,t / D / -w0.

&sDataData,@sValue
#print Started SISOARX ...
SPIN"kl;ku;1;a$,SISOARX(yd2;ud1;ml;nl;$;arx$)"
#print Stopped SISOARX ...

! Retrieve the error / residual sum-of-squares.

SSE,SCATTER(0.0;nk)
SPIN"kl;ku;1; SSE,SWAP(arx$[2+5];$)"

! Populate three parallel-arrays or data-vectors that store the (m,n,k) tuples instances.

ik,SCATTER(0;nk) 
SPIN"kl;ku;1; ik,SWAP($;$)"
&sDataData,@sValue

! Calculate the minimum SSE and its minimum index, indice, position, location, etc.

&sDataCalc,@sValue
SSEmin,MINIMUM(SSE)
SSEmini,MINIMUMI(SSE)

mi,ml
ni,nl
ki,ik[SSEmini]
&sDataCalc,@sValue
 
! Compute the "best" SISO ARX model structure found from above.
!
! * Note the arx0 data-vector is required to initially compute the w0 rate-gain for the integrating process'
!   manipulated variable, actuation or final control element.

&sDataData,@sValue
a,SISOARX(yd2;ud1;mi;ni;ki;arx0)
&sDataData,@sValue

! * FOR YOUR INTEREST * - Re-compute the a,t pseudo-white noise load disturbance using the RECURSE() data function.
!
! For interest, we can also compute the pseudo-white noise identical to thise computed by SISOARX() where we can re-compute 
! the regression residuals or prediction errors (a,t) via the RECURSE() routine below where a,t = D^2 * y,t - w0 * D * u,t-k.
!
!&sData,&sGroup
!yd2,ag
!ud1,ag
!&sData,&sGroup
!
!&sDataData,@sValue
!a,SCATTER(0.0;nt)
!a,RECURSE(ag),X1[0] - arx0[1] * X2[-ki]
!&sDataData,@sValue

! Reform the ARX model using the double differenced process output (D^2 * y,t) and the single differenced process input (D * u,t) 
! to an ARX model with no implicit differencing.

&sData,@sValue
arx, 1.0 | arx0[1] |0.0|0.0| 1 | 0 | ki |0|0|0|0|0| -99999.0 |
&sData,@sValue

! Compute the prediction of the process output y,t where y,t = yp,t + a,t or yp,t = y,t - a,t and a,t (e,t) are the regression residuals
! or prediction errors.
!
! y,t - 2 * y,t-1 + y,t-2 = w0 * u,t-k - w0 * u,t-k-1 + a,t
! y,t = 2 * y,t-1 - y,t-2 + w0 * u,t-k - w0 * u,t-k-1 + a,t
! yp,t = 2 * y,t-1 - y,t-2 + w0 * u,t-k - w0 * u,t-k-1
! a,t =  y,t - 2 * y,t-1 + y,t-2 - w0 * u,t-k + w0 * u,t-k-1

&sDataData,@sValue
yp,y-a
&sDataData,@sValue

! If the upper norm limit is zero (0.0), then do not proceed to the PID retuning.

#exit NOT(unl)

! Simulate with the actual, default or existing PID controller with the ARX dynamic model and compute
! its repsective output-error and input-move-error normalized / standardized norms (1, 2, oo).
 
&sDataData,@sValue
us,SCATTER(0.0;nt)
ys,SCATTER(0.0;nt)
&sDataData,@sValue

&sData,&sGroup
us,uys
ys,uys
&sData,&sGroup
    
!!jdk,Exxon level data "best" tuning from Kelly (1998) paper.
!&sDataData,@sValue
!pid,SWAP(-0.527;1) !Kp
!pid,SWAP( 3.425;6) !Ti
!pid,SWAP( 0.000;11)!Td
!&sDataData,@sValue
!!jdk,Unisim level data best tuning.
!&sDataData,@sValue
!pid,SWAP(      -0.300;1) !Kp
!pid,SWAP( 60.0 * 49.5;6) !Ti
!pid,SWAP(       0.000;11)!Td
!&sDataData,@sValue

&sDataData,@sValue

! * Note the degree of differencing is required below.

uys,SIMULATEPID(r;a;arx;1;pid)
&sDataData,@sValue

&sDataCalc,@sValue
es1,NORMOE(r;ys;1)
es2,NORMOE(r;ys;2)
esoo,NORMOE(r;ys;3)
dus1,NORMIME(us;1)
dus2,NORMIME(us;2)
dusoo,NORMIME(us;3)
&sDataCalc,@sValue

! Retune the PID to find the "best", smallest, lowest or minimum 1-, 2- and oo-norm output-error normalized norm subject to 
! the input-move-error norm.

&sDataData,@sValue

! * Note the degree of differencing is required below.

#print ... Started  RETUNEPID() ...
bestpid,RETUNEPID(r;a;arx;1;unl;pid)
#print ... Finished RETUNEPID() ...
&sDataData,@sValue

! Consign the total number of simulations and number that are feasible / stable i.e., no NaN, Inf nor process output 
! (y) violations or excursions.

&sCalc,@sValue
tsims,bestpid[16]
nsims,bestpid[17]
&sCalc,@sValue

&sDataData,@sValue
us2,SCATTER(0.0;nt)
ys2,SCATTER(0.0;nt)
&sDataData,@sValue

&sData,&sGroup
us2,uys2
ys2,uys2
&sData,&sGroup

! Simulate with the "best" PID settings and the SISO ARX model.

&sDataData,@sValue
pid,SWAP(bestpid[1+unt];1) !Kp
pid,SWAP(bestpid[2+unt];6) !Ti
pid,SWAP(bestpid[3+unt];11)!Td

! * Note the degree of differencing is required below.

uys2,SIMULATEPID(r;a;arx;1;pid)
&sDataData,@sValue

&sDataCalc,@sValue
ess1,NORMOE(r;ys2;1)
ess2,NORMOE(r;ys2;2)
essoo,NORMOE(r;ys2;3)
duss1,NORMIME(us2;1)
duss2,NORMIME(us2;2)
dussoo,NORMIME(us2;3)
&sDataCalc,@sValue

! Finally, add back the mean, steady-state or opening value to transform the deviation/perturbation variables back into 
! non-deviation/unperturbed variables.

&sDataData,@sValue
r,r+meanr
y,y+meany
yp,yp+meany
u,u+meanu
us,us+meanu
ys,ys+meany
us2,us2+meanu
ys2,ys2+meany
&sDataData,@sValue

! Data-groups for the "actual default", "simulated default" and simulated best" data-vector output in OML.

&sData,&sGroup
r  ,actualdefault
y  ,actualdefault
u  ,actualdefault
r  ,simulateddefault
ys ,simulateddefault
us ,simulateddefault
r  ,simulatedbest
ys2,simulatedbest
us2,simulatedbest
&sData,&sGroup