The base system consists of three areas connected in a ring configuration: Area 1 ↔ Area 2 ↔ Area 3 ↔ Area 1. Each area is represented by an aggregated model comprising a governor, a turbine, and the swing equation. Area 1 includes a thermal power plant and wind generation, along with an aggregation of electric vehicles; Area 2 consists of a hydropower plant, a PV system, and a battery energy storage system; and Area 3 comprises a gas power plant, a fuel cell, an alkaline electrolyser, and a flywheel. The synchronizing coefficients for the tie-lines are set as follows: T12 = 0.20, T23 = 0.15, and T31 = 0.18.

The Generation Dead Band (GDB) is 0.05 Hz, and the governor command limit is ±0.30 p.u. The model employs a soft limiter for the turbine power rate of change based on the TANH function; therefore, this numerical approximation should be distinguished from a hard limiter. All states are initialized to zero at t = 0. An open-loop simulation was conducted for 60 seconds on the 6001-bus network, with an output step size of 0.01 seconds. The ode45 solver was used for the delay-free case, while a delay-aware solver was employed for the case involving delays; a standalone Simulink implementation—utilizing ode45 and delay blocks—was also evaluated.

The state vector for each area includes frequency deviation, turbine mechanical power, and the governor command:

x = [Δf , ΔP , ΔP ] , i ∈ {1,2,3}. (1)

i i m,i gov,i

The swing equation incorporates renewable power generation, load, interconnector exchange, and damping:

dΔfi

In the step-load scenario, the maximum absolute frequency deviations in Areas 1, 2, and 3 were 0.1815888696 Hz, 0.1607677882 Hz, and 0.1687669125 Hz, respectively, while the maximum tie-line power deviation was 0.05400969053p.u. At [time to be specified], the frequency deviations in Areas 1, 2, and 3 were −0.1275122359 Hz, −0.1108743745 Hz, and −0.1313356948 Hz, respectively. The final discrepancy was 0.0204613204 Hz. Because secondary Automatic Generation Control (AGC) is inactive in the baseline model, the system state at [time to be specified] cannot be regarded as a verified steady state.

In the wind scenario, the maximum frequency deviations in the three areas were 0.1207079814Hz, 0.1309725371Hz, and 0.1358342313Hz, and the maximum tie-line power was reported as 0.03898786106p.u. The final values were -0.03388383831Hz, -0.07926832381Hz, and -0.08834699974Hz, respectively, with a final discrepancy of 0.0544631614Hz. The oscillatory and bounded nature of the wind disturbance results in a response that is smoother yet more oscillatory. In the PV step scenario, the maximum absolute frequency deviations in areas 1 through 3 were 0.1397334902, 0.1905673228, and 0.1473227631 Hz, respectively, while the maximum tie-line power was 0.06379595528p.u. The final values were +0.1185087119, +0.06499525771, and +0.1313140045Hz, with a final deviation of 0.0663187468 Hz. The positive sign is consistent with the interpretation that an increase in PV generation acts as a net reduction in load.

In the combined scenario, the maximum frequency deviations were 0.1852248572, 0.2038793548, and 0.2167180939Hz, and the maximum tie-line power was 0.1063107815p.u. The final values were - 0.06310970256, -0.06719563367, and -0.1036986276 Hz, with a final deviation of 0.040588925 Hz. Consequently, the combined scenario produced the highest single-area peak for the base model as well as the greatest stress on the tie-line.

In the single-area model with a load of 0.1 p.u., the analytical value for the steady-state frequency—accounting for the Governor Dead Band (GDB)—was -0.279922779923 Hz. ode45 yielded a value of -0.279922779984 Hz, with an absolute error of 6.13663009297×10⁻¹¹ Hz—a value significantly smaller than the 10⁻³Hz limit. In the three-area model, the analytical equilibrium was -0.124629196 Hz, while the numerical average at t=60 s was -0.123240768 Hz; the difference is 0.00138842712 Hz, which exceeds the 10⁻³ Hz agreement threshold. This discrepancy must be explicitly acknowledged and should not be presented as a perfect match.

In the single-area model with a 0.1 p.u. load, the analytical steady-state frequency—accounting for GDB— was -0.279922779923 Hz; `ode45` yielded a value of -0.279922779984 Hz, with an absolute error of 6.13663009297×10⁻¹¹ Hz—a value significantly smaller than the 10⁻³ Hz limit. In the three-area model, the analytical equilibrium was -0.124629196 Hz, while the numerical average at t=60s was -0.123240768 Hz; the difference is 0.00138842712 Hz, which exceeds the 10⁻³ Hz agreement threshold. This discrepancy must be explicitly acknowledged and should not be presented as a perfect match.

In an independent comparison of Simulink, the maximum absolute difference relative to MATLAB was

−7 −5

reported as 2.4 ∗ 10 for the load-step scenario and 9.2 ∗ 10 for the delay sweep. These findings support the numerical consistency between the two implementations within the specified tolerances; however, the discrepancy regarding the terminal frequency criterion in the three-area model persists.

For Scenario 1, a delay of τ ∈ 0.50, 100, 150, 200 ms was applied. The maximum Δf1 increased as the delay increased, rising from 0.18158887 to 0.20567723, 0.50758020, 0.89603921, and 1.10980591 Hz; the maximum values across all areas were 0.18158887, 0.35942208, 0.67855475, 1.32501995, and 1.62017753 Hz, respectively. At larger delays, a sign change in the final frequency value indicates oscillatory behaviour induced by the delay. These results demonstrate the necessity of explicitly accounting for delay when evaluating the controller.

The fuzzy input vector consists of the control error E and the change in error CE. Five linguistic labels— NB, NS, Z, PS, and PB—are defined, with centers at -1.00, -0.50, 0.00, 0.50, and 1.00; a width of w = 0.30; and a mean uncertainty half-width of Δm = 0.10. Four α levels are set at 0.25, 0.50, 0.75, and 1.00, and the MacVicar–Whelan rule base comprises 25 rules. The maximum value of the UMF is reported as 1.000, while the maximum value of the LMF is 0.667.

|x − (c − Δm)| |x − (c + Δm)|

t1(x) = max(1 − , 0) , t2(x) = max(1 − , 0). (7)

UMF(x) = max(t , t ), LMF(x) = min(t , t ), H(x) = [UMF(x) + LMF(x)]. (8)

1 2 1 2

loα(x) = (1 − α)LMF(x) + αH(x), upα(x) = (1 − α)UMF(x) + αH(x). (9)

∑a αa yα

y = . (10)

∑a αa

In each α-plane, the lower and upper membership degrees are calculated as a combination of the LMF/UMF and an intermediate function H; the final output is obtained via a weighted average based on the α-level. Due to uncertainty in the mean, the UMF of each label may exhibit two symmetric peaks flanking the nominal centre; this characteristic is preserved within the membership function data.

1 1 1 2 3 1 1 2 3 4

rulemap = 1 2 3 4 5 . (11)

2 3 4 5 5 [3 4 5 5 5]

yconseq = {−1.0, −0.5,0.0, +0.5, +1.0}, fr,l = min(loα(E), loα(CE)), fr,u =

(12)

min(upα(E), upα(CE)).

At the point E=0.30, CE=0.10, and α=0, only two rules—with the values shown in the figure—are activated: Rule 13, with a lower strength of 0.3333 and an upper strength of 0.6667; and Rule 18, with a lower strength of 0 and an upper strength of 0.6667. In the reference test case where α=0.5, the strengths are as follows: f13,l=1/12, f13,u=1/4, f18,l=1/6, and f18,u=1/2.

In the NT method, the output for each α-level is the weighted average of the rule consequences, using the sum of the lower and upper strengths. If the denominator is less than 10⁻¹², an output of zero is returned.

∑r=1(fr,l + fr,u) yconseq,r

NT

yα = . (13)

∑ (fr,l + fr,u)

In the KM method, the lower and upper bounds are calculated by iterating over the switching point k within

yl+yu

the sorted consequences. Subsequently, the midpoint is taken as the output for a given α-level, and the α-weighted average is used as the final output.

∑r≤k fr,u yr + ∑r>k fr,l yr

yl = . (14)

∑r≤k fr,u + ∑r>k fr,l

∑r≤k fr,l yr + ∑r>k fr,u yr

yu = . (15)

∑r≤k fr,l + ∑r>k fr,u The fractional block FR is applied to the output time series NT. The weight of each sample is defined by the time interval and the fractional order β, while an exponential factor with a constant Tf = 0.1s limits the memory:

β−1 tk − tj ∑j<k wj yNT(tj)

w = (t − t ) exp [− ] , y (t ) = . (16)

j k j mem k

Tf ∑j<k wj

yFR(tk) = (1 − β)ymem(tk) + βyNT(tk). (17)

At β=1, the FR block reduces exactly to NT; as β approaches zero, the contribution of the causal memory increases. The selection policy is as follows: The NT method for τd < 0.05s; the FR method with β=0.5 for 0.05 ≤ τd < 0.15s; and the KM method for τd ≥ 0.15s. While this rule dictates the algorithm selection, it does not, by itself, guarantee closed-loop stability or robustness against large delays.

NT, τd < 0.05 s,

select(τd) = {FR (β = 0.5), 0.05 ≤ τd < 0.15 s, (18)

KM, τd ≥ 0.15 s.

Three independent tests were conducted to verify the implementation of the control method. First, the maximum absolute error between the FR method (with β = 1) and the NT output was found to be 0.000000*100. Second, at the point E, CE, α = 0.30,0.10,0.50, both the direct counting of switching points and the iterative KM algorithm yielded yl = 0.20000000 and yu = 0.42857143 = 3/7; the difference between the two bounds was zero. Third, due to the symmetry of the rule base, the outputs of both the NT and KM methods at the origin E, CE = 0,0 were +0.00000000.

In a static sweep with CE = 0 and 201 points ranging from -1 to +1, the mean value of yNT−E was 0.0527, while the mean value of yKM−E was 0.0497. The maximum difference between the NT and KM outputs was 0.0097, and the output bounds for both methods were 1.0000pm. A dynamic test was performed using E(t) = 0.8sin2π ⋅ 0.3t and CE(t) = [Et(k) − Et(k − 1)]/Δt, with Δt = 0.05s over the interval [0, 10]s. The Root Mean Square Error (RMSE) was 0.0045 between the NT and KM methods, 0.2527 between the NT and FR methods (with β =0.5), and 0.3438 between the NT and FR methods (with β = 0.2). This discrepancy reflects the effect of fractional memory, rather than necessarily an improvement in the final control criterion. The reported computational time for each evaluation is 14.987792 ms for NT and 15.247681 ms for KM, representing a 1.02 ratio. These figures are specific to the implementation and hardware tested; they should not be generalized to other systems without information regarding an identical computing environment.

The three raw Pareto archive objectives are Integral of Time-weighted Absolute Error (ITAE), Integral of Squared Error (ISE), and Total Variation (TV). ITSE is also used for the purpose of a comparative report:

ITAE = Δt ∑ t [∑|Δf (k)| + ∑|ΔP (k)|]. (19)

k i tie,i k i=1 i=1

2 2 (20)

ISE = Δt ∑ [∑ Δ fi(k) + ∑ Δ Ptie,i(k) ].

k i=1 i=1 2 (21)

ITSE = Δt ∑ tk ∑ Δ fi(k) , TV = ∑ ∑ |ui(k + 1) − ui(k)|.

k i=1 i=1 k

The strict constraints relate to the frequency deviation |Δf| ≤ 0.20 Hz and the tie-line power deviation |ΔPtie| ≤ 0.10 p.u. An incremental penalty with a coefficient of 105 was applied:

Ppen = 10 [max(0, max|Δf| − 0.20) + max(0, max|ΔPtie| − 0.10)]. (22)

F = [ITAE + Ppen, ISE + Ppen, TV + Ppen]. (23)

In GA and PSO, the weighted cost is formulated based on three raw objectives plus the aforementioned penalty term. The uncontrolled baseline case was infeasible due to severe constraint violations; it recorded a cost of 895,105.929999, a maximum Δf of 6.505847 Hz, and a maximum ΔPtie of 2.741728 p.u.

The MOPSO archive contained 7 solutions. The ranges were reported as follows: ITAE [80.2939, 83.8939], ISE [0.65561, 0.75757], and TV [23.7056, 26.5612]. A compromise solution—the 5th out of the 7—was selected, having a normalized L2-distance of 0.7145 from the ideal point. Its parameter vector is: Θ*=0.000000, 1.871288, 2.000000, 0.000000, 1.213327, 1.000000, 0.679925, 0.571083.

The corresponding values were: ITAE = 81.722190, ISE = 0.6906788, TV = 25.089462, maximum |Δf| = 0.169317 Hz, maximum |ΔPtie| = 0.053679 p.u., and the feasibility flag was 1.

In the weighted-cost comparison, GA yielded ITAE = 41.286386, ISE = 0.2785581, ITSE = 2.494403, and TV = 45.969514; the maximum values were |Δf| = 0.197900 Hz and |ΔPtie| = 0.075064 p.u. PSO yielded ITAE = 34.085742, ISE = 0.1708558, ITSE = 1.655889, and TV = 42.675775; the maximum values were Δf = 0.145928 Hz and |ΔPtie| = 0.062967 p.u. Both solutions were feasible within the reported constraints. These values alone do not establish a definitive, unified ranking among MOPSO, GA, and PSO, as MOPSO optimized a multi-objective vector, whereas GA and PSO optimized a weighted-cost function; this fundamental difference in the objective function must be taken into account when interpreting the comparison.

Over 25 iterations, the best ITAE reported for MOPSO decreased from 265.967268 to 80.293942; the final archive size was 7, representing a relative reduction of 69.81%. The cost for GA decreased from 27.615303 to 25.819880 (6.5%), while the cost for PSO decreased from 29.692369 to 22.237794 (25.1%). In 10 independent runs of PSO, the following statistics were recorded: a mean of 31.1937, a standard deviation of 6.3186, a minimum of 22.9568, a maximum of 37.9097, and a coefficient of variation of 20.26%. Consequently, the results exhibit significant dispersion, and the outcome of a single run should not be presented as the definitive result for all runs.

According to Fig. 21, the range for ITAE is [80.2939, 83.8939], for ISE is [0.65561, 0.75757], and for TV is [23.7056, 26.5612]; a compromise solution is indicated.

According to Fig. 22, the correlation between ITAE and ISE is reported as 0.997, the correlation between ITAE and TV as -0.938, and the correlation between ISE and TV as -0.944.

According to Fig. 26, the peak values of Δf1 for MOPSO, GA, PSO, and the baseline are 0.1605, 0.1778, 0.1251, and 3.3195, respectively. The final values of Δf1 are -0.10659, -0.02544, -0.00983, and -0.68886, respectively.

According to Fig. 27, the mean is 31.1937, the standard deviation is 6.3186, and the coefficient of variation is reported as 20.26%.

The execution-data archives indicate that fractional-order memory affects the closed-loop response, as the ITAE values obtained with the FOGT2-Adp controller differ from those of FOGT2-b1 in all four scenarios. Specifically, the differences, defined as ITAE(FOGT2-Adp) − ITAE(FOGT2-b1), are [0.71483701, 2.4122148, 2.8655018, 3.429844]. These positive values indicate that FOGT2-Adp yields a higher ITAE than FOGT2-b1 in every scenario under the current calibration. Furthermore, the differences ITAE(PID) − ITAE(FOGT2-Adp) are [−5.7944549, 8.2061213, −3.9289591, 9.5366412]. Thus, the PID controller achieves a lower ITAE than FOGT2-Adp in Scenarios S1 and S3, whereas FOGT2-Adp achieves a lower ITAE in Scenarios S2 and S4. These results do not support the claim that the proposed controller provides the best performance across all scenarios. Instead, they indicate that the relative performance of the controllers depends on the disturbance scenario and the evaluation metric.

In the selected calibration, the PID controller violated both defined limits, with a maximum |Δf| of 0.2366023 Hz and a maximum |ΔPtie| of 0.10247591. FOGT2-Adp also exceeded the 0.20 Hz limit in scenario S3, reaching a maximum |Δf| of 0.20126397 Hz. These values should not be truncated or omitted; the violation of constraints constitutes a valid part of the scientific findings. Furthermore, the results for FOGT2-b1 and IT2-FPID in the primary tables are identical, which is consistent with setting β=1 and applying an α-plane of 0 in the reported implementation. This structural equivalence should not be interpreted as two distinct, independent results.

In the data archive for τ = 150 ms and τ = 200 ms, all five controllers were deemed invalid, exhibiting infinite ITAE values and NaN (Not a Number) steady-state values; consequently, the data do not support a stable delay tolerance of up to 200 ms. At τ = 100 ms, frequency peaks in the range of 103 to 104 Hz were observed. The adaptive type-reduction selection strategy, on its own, failed to resolve this performance breakdown. In the Monte Carlo simulation, 100 out of 100 runs for the PID controller and 100 out of 100 runs for the FOGT2-Adp controller were recorded as valid. The mean ITAE was 126.29831 (with a standard deviation of 1.2431464) for the PID controller, and 131.45981 (with a standard deviation of 1.2620647) for the FOGT2- Adp controller. In this specific test set, the comparison favours the PID controller in terms of mean ITAE; however, the number of valid runs alone does not necessarily imply superior performance quality.

In the linearized closed-loop pole analysis, the PID controller was found to be unstable, with the maximum real part of its poles equal to 0.005813647. In contrast, the T1-FPID, IT2-FPID, FOGT2-b1, and FOGT2- Adp controllers satisfied the stability criterion applied to the linearized model. For the FOGT2-Adp controller, the Gain Margin (GM) was [value to be specified], while the phase margin (PM) was 64.228196°. These findings should be interpreted with caution, as the linearized analysis does not fully capture the effects of nonlinearities and dynamic characteristics, including saturation, dead areas, rate limits, fractional-order memory, and communication delays. Consequently, the results do not establish the overall stability of the nonlinear closed-loop system under time delays and actuator constraints; this requires further analysis using methods that explicitly account for these effects.

Performance metrics reported in related studies also vary considerably with respect to system configuration, modeling assumptions, disturbance conditions, controller architecture, and evaluation criteria. For example, the results reported for islanded microgrid controllers [1], hybrid power systems [2], parallel control architectures [3], adaptive fractional-order controllers [4], Type-2 fuzzy controllers for nonlinear systems [9], general Type-2 fractional-order PID controllers [5], interconnected microgrids [8], systems subject to time delays [6], and adaptive Type-2 control strategies [7] were obtained under different experimental or simulation conditions. Accordingly, the supplementary comparison table is intended to contextualize the proposed framework in relation to the architectures and methodological approaches adopted in previous studies, rather than to establish a direct quantitative ranking of controller performance. Similarly, Darvish Falehi [10] investigated a fractional-order Type-2 fuzzy control structure for LFC using multi-objective optimization. However, the numerical results of that study cannot be directly compared with those obtained in the present work unless the underlying models, objective functions, constraints, disturbance scenarios, and evaluation protocols are sufficiently aligned.