Optimal tuning of multi-PID controller using improved CMOCSO algorithm

Optimization objectives

PID controllers, which are designed to minimize time-domain performance metrics, have been extensively employed in control engineering. Historically, some researchers have concentrated on optimizing PID parameters with a singular objective; however, this methodology is inadequate for identifying the optimal PID parameters. Therefore, a holistic assessment of multiple factors is essential. This article considers several standard performance metrics simultaneously to evaluate the efficacy of electro hydraulic servo synchronous control systems.

The system overshoot quantifies the percentage by which the peak value of the step response surpasses its final steady-state value. It is computed as illustrated in Eq. (3).

(3) $${\mathrm{J}}_1=\sigma =\sum _{i=1}^n\frac{c_i^{max}-c_i(\mathrm{\infty })}{c_i(\mathrm{\infty })},$$where $c$ represents the step response curve of the system, $c_{max}$ denotes the maximum value of the step response, $c(\mathrm{\infty })$ denotes the final steady state value of the step response, $c_i$ denotes the step response curve associated with the i-th PID controller within the system, and $n$ signifies the number of PID controller subsystems in the system.

The adjustment time is defined as the minimum duration required for the step response to attain and remain consistently within a 2% error margin of the final steady-state value. The calculation is shown in Eq. (4).

(4) $${\mathrm{J}}_2=t_s=\sum _{i=1}^n\left\{t|\left|\frac{c_i(t)-c_i(\mathrm{\infty })}{c_i(\mathrm{\infty })}\right|<0.02\right\},$$where $c$ represents the step response curve of the system, $c(\mathrm{\infty })$ denotes the final steady state value of the step response, $c_i$ denotes the step response curve associated with the i-th PID controller within the system, $t$ denotes the system’s simulation time, and $n$ stands for the number of PID controller subsystems.

The rise time is defined as the duration necessary for the step response to first reach the final steady-state value. It is calculated as shown in Eq. (5).

(5) $${\mathrm{J}}_3=t_r=\sum _{i=1}^n\left\{t|c_i(t)-c_i(\mathrm{\infty })=0\right\},$$where $c(\mathrm{\infty })$ denotes the final steady state value of the step response, $c_i$ denotes the step response curve associated with the i-th PID controller within the system, $t$ denotes the simulation time of the system, and $n$ denotes the number of PID controller subsystems in the system.

The peak time is defined as the time required for the step response to surpass the final steady-state value for the first time. It is calculated as shown in Eq. (6).

(6) $${\mathrm{J}}_4=t_m=\sum _{i=1}^n\left\{t|c_i(t)>c_i(\mathrm{\infty })\right\},$$where $c(\mathrm{\infty })$ denotes the final steady state value of the step response, $c_i$ ddenotes the step response curve associated with the i-th PID controller within the system, $t$ denotes the simulation time of the system, and $n$ denotes the number of PID controller subsystems in the system.

The steady-state error quantifies the difference between the theoretical and actual step response values, computed as outlined in Eq. (7).

(7) $${\mathrm{J}}_5=se=\sum _{i=1}^n\left|out-c_i(\mathrm{\infty })\right|,$$where $c(\mathrm{\infty })$ denotes the final steady state value of the step response, $c_i$ denotes the step response curve associated with the i-th PID controller within the system, $out$ is the ideal output value of the system, $t$ denotes the simulation time of the system, and $n$ denotes the number of PID controller subsystems in the system.

The synchronization error assesses the deviation between the master and slave structures in systems utilizing master-slave synchronization control, calculated according to Eq. (8).

(8) $${\mathrm{J}}_6=e_{ss}=\sum _{i=1}^{n-1}\int c^m(\mathrm{\infty })-c_i^s(\mathrm{\infty })dt,$$where $c(\mathrm{\infty })$ denotes the final steady state value of the step response, $c_i$ denotes the step response curve associated with the i-th PID controller within the system, $t$ denotes the simulation time of the system, $n$ denotes the number of PID controller subsystems in the system, $c^m$ denotes the step response curve of the PID controller of the master structure, and $c^s$ denotes the step response curves of the PID controllers of the slave structure.

The integral time squared error not only mitigates excessively large signal values that may exceed the amplitude of the PID controller, but also avoids the error rate that may cause delays in the sensors in the controller (Lin et al., 2022). The calculation is illustrated in Eq. (9).

(9) $${\mathrm{J}}_7=\mathrm{I}\mathrm{T}\mathrm{S}\mathrm{E}=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\int }_0^{\mathrm{\infty }}t{e}_i^2(t)dt,$$where $e(t)$ denotes the error between the ideal and actual output values of the system at time $t$, $i$ denotes the i-th PID controller subsystem within the system, $t$ signifies the simulation time, and $n$ denotes the number of PID controller subsystems in the system.

Constraints set

In order to define parameters for system robustness, two essential factors are taken into account: amplitude margin and phase margin. The phase margin (PM) refers to the permissible additional phase lag at the gain crossover frequency that is necessary to achieve critical stability within the system. Conversely, the gain margin (GM) represents the maximum allowable increase in the system’s gain while still preserving critical stability. As outlined in Hu (2013), the PM is generally observed to fall within the range of 30 to 60 degrees, whereas GM is typically expected to exceed 6 dB.

The magnitude of system overshoot is a critical consideration in the assessment of system stability. An excessive overshoot may result in instability and adversely affect the system’s normal functioning. Consequently, it is imperative to establish a constraint on the level of overshoot. By regulating overshoot, one can enhance system stability during operation and mitigate unwarranted fluctuations and atypical behaviors.

The system overshoot constraint calculation is calculated by Eq. (10).

(10) $${\sigma }^{\mathrm{c}}=\underset{i=\mathrm{1..}n}{max}\left(\frac{c_i^{max}-c_i(\mathrm{\infty })}{c_i(\mathrm{\infty })}\right),$$where $c$ denotes the step response curve of the system, $c^{max}$ denotes the maximum value of the step response, $c(\mathrm{\infty })$ denotes the final steady-state value of the step response, $c_i$ denotes the step response curve associated with the i-th PID controller within the system, and $n$ denotes the number of PID controller subsystems in the system.

Individual coding and parameter search range

The ‘Multi-PID controller model’ discussed in the preceding section illustrates that the simulation framework comprises four PID controllers, each characterized by three parameters: proportional, integral, and derivative. Consequently, the decision space for the constrained multi-objective optimization is represented in 12 dimensions. Specifically, the decision variables are organized as [ $K_{p1}$, $K_{i1}$, $K_{d1}$, $K_{p2}$, $K_{i2}$, $K_{d2}$, $K_{p3}$, $K_{i3}$, $K_{d3}$, $K_{p4}$, $K_{i4}$, $K_{d4}$], where $K_{pj}$, $K_{ij}$, and $K_{dj}$ denote the proportional, integral, and derivative parameters of the j-th PID controller, respectively. The decision space is encoded using real numbers, with each PID parameter expressed as a real number. The search space is established based on parameters referenced in the literature (Zhang et al., 2016) and is defined by extending the range both upwards and downwards by a specified multiple (for instance, from 0.5 to 1.5 times the original value). For instance, if the PID parameter value calculated from the literature is denoted as X, the search range can be delineated as $[0.5X,1.5X]$. Should it be observed that the majority of the optimized PID parameter results tend to cluster near the boundaries of the search range, the multiplicative factor of the search range may be further modified to broaden the search space.

In summary, the constrained multi-objective optimization model proposed in this article is calculated by Eq. (11).

(11) $$\begin{array}{c}minf(x)=({\mathrm{J}}_1,{\mathrm{J}}_2,{\mathrm{J}}_3,{\mathrm{J}}_4,{\mathrm{J}}_5,{\mathrm{J}}_6,{\mathrm{J}}_7),\hfill \\ x=[K_{p1},K_{i1},K_{d1},K_{p2},K_{i2},K_{d2},K_{p3},K_{i3},K_{d3},K_{p4},K_{i4},K_{d4}],\hfill \\ s.t.GM>g_{min},{\mathrm{p}}_{max}>PM>p_{min},{\sigma }_{max}>{\sigma }^{\mathrm{c}},\hfill \end{array}$$where $J_1$ is the overshoot, ${\sigma }^c$ is the overshoot constraint, ${\sigma }_{max}$ is the upper bound of the overshoot, $J_2$ is the rise time, $J_3$ is the regulation time, $J_4$ is the peak time, $J_5$ is the steady state error, $J_6$ is the error between the two subsystems, $J_7$ is the integral squared error, GM and PM denote the gain margin and phase margin of the system, respectively, and $g_{min}$, $p_{max}$, $p_{min}$ denote the upper and lower bounds of the respective constraints.

Two-stage policy

ICMOCSO is an improved CMOCSO algorithm that utilizes a two-stage strategy. In the traditional CMOCSO, two populations coexist throughout the evolutionary process, resulting in increased algorithm runtime. However, ICMOCSO divides the entire evolutionary process into two stages by introducing a centroid movement strategy. In the initial phase, a new grouping strategy and environment selection method are implemented, accelerating the convergence of the main_Pop. Simultaneously, coevolution is achieved by sharing the offspring of main_pop and help_pop. In the second phase, emphasis on the evolution of the main_Pop further accelerates the entire algorithm. The primary principle of the ICMOCSO algorithm is illustrated in Fig. 3.

Combined two-stage strategy

In traditional CMOCSO, the presence of two continuously evolving populations throughout the evolutionary process often results in prolonged runtime and increased computational costs when addressing real-world problems. Conversely, ICMOCSO algorithm is structured into two distinct phases. The initial phase introduces a variety of novel strategies, including grouping strategies and environmental selection methods, aimed at enhancing the convergence of the main_Pop while facilitating coevolution with the help_Pop. During this phase, the main_Pop employs $\epsilon$-constraint techniques (Song et al., 2024) to conduct a global search for the CPF, while the help_Pop focuses on identifying the PF within the objective space. This strategic design mitigates the risk of the main_Pop prematurely converging to local optima, thereby promoting the global convergence of the algorithm. In the environmental selection process, the main_Pop selects optimal solutions from the offspring generated by both the main_Pop and help_Pop, ensuring the continued evolution of the main_Pop. Once the centroid of the main_Pop reaches a predetermined distance, the algorithm transitions into the second phase.

In the second phase, since the first phase has already determined the approximate range of the CPF, only the main_Pop continues to evolve under all constraints. The focus of this phase is to accelerate convergence and reduce the algorithm’s cost. By focusing on the evolution of the main_Pop, the algorithm can find the optimal solutions faster, effectively shortening the runtime while ensuring search accuracy and effectiveness.

The $ϵ$-constraint technique is an approach utilized to manage constraints dynamically by incorporating a relaxation factor, $ϵ$, within the optimization process. This methodology permits the initial relaxation of stringent constraint requirements during the early phases of the optimization search, thereby broadening the search space and augmenting the diversity of potential solutions. As the optimization process advances, the value of $ϵ$ is systematically reduced, which progressively strengthens compliance with the constraints. This mechanism effectively directs the algorithm towards convergence within the feasible solution region that fully adheres to the established constraints. As outlined in reference Fan et al. (2019), the updating process for $ϵ$ is as follows:

(12) $$\begin{array}{c}ϵ(t)=\{\begin{array}{ccc}{\varphi }^t(x),\hfill & if& t=0,\hfill \\ (1-\tau )ϵ(t-1),\hfill & if& r_t<\alpha ,\hfill \\ {\varphi }_{max}{(1-\frac{t}{T_c})}^{cp},\hfill & if& r_t,\ge \alpha ,\hfill \\ 0,\hfill & if& t\ge T_c,\hfill \end{array}\hfill \end{array}$$where ${\varphi }^t(x)$ is the total constraint violation of the top individuals in the population, $r_t$ is the feasible ratio of population, $\alpha$ and $\tau$ are two control parameters, ${\varphi }_{max}$ is the maximal constraint violation of population found so far, and $T_c$ controls the maximal generation that $ϵ$ works, while $cp$ reducing the $ϵ$ value in the case of $r_t\ge \alpha$.

Center-point movement strategy

Currently, many two-phase algorithms use a semi-partitioning method to distinguish stages, considering the first $N\mathrm{\%}$ of the entire search process as one phase and the remaining part as another phase. However, the ICMOCSO algorithm differs from this common semi-partitioning method in that it features an automatic phase transition mechanism. This mechanism determines the current phase by monitoring the movement distance of the central point of the main population. When the movement distance of the main population’s central point is sufficiently small, it indicates that the population has found a relatively stable optimal solution in the solution space and the evolutionary process has reached a saturation state. At this point, continuing the search in the current phase may not yield further improvements, so the algorithm can automatically transition to the next phase to better explore other parts of the solution space. The central point of the population aggregates and summarizes the information of all individuals in the population. By analyzing the position and movement of the central point, one can obtain global information about the population in the solution space. The calculation formula for the central point is calculated by Eq. (13).

(13) $$C_t=\frac{1}{\left|MainPop\right|}\sum _{x\in MainPop}F(x),$$where C represents the central point of the main population, $t$ represents the number of iterations of the algorithm, $|MainPop|$ represents the number of individuals in the main population, $x$ represents an individual in the main population, and $F(x)$ represents the fitness value of each individual.

If the Euclidean distance between the central points of the main population in two iterations exceeds a critical value, the search process will remain in the first phase; otherwise, it will automatically transition to the next phase. The formula for the central point movement strategy is calculated by Eq. (13).

(14) $$\{\begin{array}{c}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}1:\left||{\mathrm{C}}_t-C_{t-1}\right||>\theta ,\hfill \\ \mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}2:\left||{\mathrm{C}}_t-C_{t-1}\right||\le \theta ,\end{array}$$where C represents the central point of the main population, $t$ represents the number of iterations of the algorithm, $\theta$ represents the minimum value of the central point movement, $stage1$ and $stage2$ represent the two phases of the improved CMOCSO, and || || calculates the Euclidean distance between the two central points.

Grouping strategy

The ICMOCSO algorithm presents an innovative grouping strategy that integrates pareto dominance with fitness value comparison. In contrast to the CMOCSO approach, where individuals are randomly grouped and subsequently paired for mutual learning based on their fitness values, the ICMOCSO method categorizes particles into winning and losing groups through the application of pareto dominance and fitness value assessments. In this framework, non-dominated individuals are immediately classified as winners. Among the remaining individuals, those exhibiting lower fitness values are selected to join the winning group, thereby maintaining an equal distribution of individuals across both groups. Following the completion of this comparison, the winners advance directly to the subsequent generation, while the losers engage in a learning process from the winners. This strategic approach not only improves the search efficiency of the algorithm but also enhances the overall quality of the solutions produced. The individual comparison process is shown in Fig. 4.

The steps of the grouping strategy in this algorithm are as follows: Firstly, individuals in the population are divided into different front groups through non-dominated sorting. Then, individuals are classified into winner or loser groups based on whether they are in the first layer of the pareto front. If no individuals are in the first layer of the pareto front, individuals are allocated to the two groups based on fitness value comparison. If the number of individuals on the pareto front exceeds half of the total population, the top $N/2$ individuals with the smallest fitness values are selected as the winner group, and the remaining individuals are the loser group. If the number of individuals on the pareto front is less than half of the total population, the top $N/2$ individuals with the highest fitness values are selected as the loser group, and the remaining individuals are the winner group.

Each individual is allocated a fitness value, which is determined by two primary factors: the amount of the solution dominated by other solution $ϵ$ constraints, and the crowding distance that measures the proximity of the solution to other solutions. The precise methodology for this calculation is presented in Eq. (15).

(15) $$F(i)=|\{j\in Q|j{\prec }_{ϵ}i\}|+d,$$where || denote the cardinality of a set, $d$ signify the crowding distance between two solutions, Q represent the current solution set, and the notation $j{\prec }_{ϵ}i$ indicate that solution $j$ $ϵ$-constraint dominates solution $i$. For two solutions, $x$ and $y$, it is stated that $x$ $ϵ$-constraint dominates solution $y$ (denoted as $x{\prec }_{ϵ}y$) if one of the following conditions is satisfied (Takahama, 2006):

• $\varphi (x)<ϵ\text{\&}\varphi (y)<ϵandx\prec y$,

• $\varphi (x)=\varphi (y)=ϵandx\prec y$,

• $\varphi (x)>ϵ\text{\&}\varphi (y)>ϵand\varphi (x)\prec \varphi (y)$,

where $ϵ$ denotes relaxation factor and $\varphi (x)$ denotes overall constraint violation of $x$, $x\prec y$ denotes that $x$ pareto dominates $y$. $x\prec y$ means that $x$ is not worse than $y$ on all objective functions, and $x$ is better than $y$ on at least one objective function.

Environmental selection strategies

The conventional environmental selection strategy employed in CMOCSO relies on the Euclidean distance archive truncation method. While the Euclidean distance metric offers certain benefits in quantifying spatial relationships among individuals, it fails to account for the distribution characteristics of the solution space. In particular, the distribution of individuals within the solution space may be uneven, a nuance that the euclidean distance does not adequately address. In contrast, the shifted density estimation (SDE) strategy (as referenced in Li, Yang & Liu, 2013) provides a more comprehensive evaluation of solution quality by considering both convergence and diversity. This approach has gained considerable traction in various multi-objective evolutionary algorithms (MOEAs) (as noted in Bai & Zhang (2024), Ou et al. (2024)). The fundamental premise of the SDE strategy is to assess the quality of an individual within the population based on the density of neighboring individuals within a specified spatial range. Specifically, the SDE strategy emphasizes the density surrounding an individual; a high surrounding density suggests that the solution represented by that individual is situated in a relatively crowded region, potentially indicating its superiority. Conversely, a low surrounding density may imply that the solution is located in a relatively sparse area, signaling a need for further exploration or enhancement.

The environmental selection process of the ICMOCSO algorithm follows these steps: First, calculate the fitness values of all individuals. Individuals with fitness values less than 1 are selected as the next generation population. If the number of the next generation population is equal to the preset population size, the environmental selection process ends. If the next generation population size is less than the preset population size, individuals with smaller fitness values are selected from unselected individuals to join the next generation population. If the next-generation population size exceeds the preset population size, individuals with higher shifted density are deleted until the next generation population size reaches the preset size.

Complexity analysis

As previously discussed, the proposed ICMOCSO framework encompasses several key components, including the central-point moving strategy, grouping strategy, $stage1$, $stage2$, and the environment selection strategy. The time complexity associated with the central-point moving strategy is primarily influenced by the number of fitness evaluations, center-point calculations, and stage transition assessments, resulting in a complexity of $O(N)$. In contrast, the grouping strategy exhibits a time complexity of $O(N^2)$, which is largely determined by non-dominated sorting, as well as fitness value comparisons and calculations. The computational demands of the CMOCSO utilized in $stage1$ are characterized by a complexity of $O(N^3)$, while the main_Pop update employed in $stage2$ incurs a complexity of $O(m{N}^2)$, as detailed in the original literature. Regarding the environment selection strategy, its fitness calculations and truncation operations are analogous to those found in SPEA2, thereby yielding a computational complexity of $O(N^2)$. In summary, the overall computational complexity of the ICMOCSO framework is $O(N^3)$.

Proof of algorithm effectiveness

The experiments were conducted on a PC with Windows 10 operating system, an Intel(R) Core(TM) i5-9500 CPU @ 3.00 GHz, and 8 GB RAM, using MATLAB 2021b. The multi-objective optimization algorithm testing platform, PlatEMO (Tian et al., 2017) (version 4.2), was employed for performance comparisons.

To evaluate the performance of the ICMOCSO algorithm, we compared it with several other algorithms: CTSEA (Ming et al., 2021), CMOCSO (Ming et al., 2022a), MaOEAIT (Sun et al., 2019), RVEAa (Cheng et al., 2016), TriP (Ming et al., 2022b), and TiGE2 (Zhou et al., 2018). Each algorithm was run 30 times independently on 16 ZXH_CF (Zhou, Xiang & He, 2020) test sets to ensure the stability and reliability of the results. For fair performance comparisons, we used the hyper volume (HV) metric (Riquelme, Von Lücken & Baran, 2015) as the evaluation standard. HV is widely used in the field of multi-objective optimization as it comprehensively reflects the distribution and uniformity of the solution set. In the conducted experiment, we established several public parameters: the maximum number of iterations was set to 10,000, the population size was 100, the number of targets was determined to be 7, and the number of decision variables was set at 12. Additionally, the proposed algorithm was configured with a minimum value for the two-stage switching center shift of 0.005, a rate of decrease for the $ϵ$ value ( $cp$) was controlled at a value of 2, a maximal generation ( $T_c$) of 900, and values for $\alpha$ and $\tau$ set at 0.95 and 0.05, respectively. For other algorithms, the parameters are set according to the original article, and the settings of public parameters are consistent with those in this article.

The HV metric results for the ZXH_CF test functions are detailed in Table 2. In Table 2, “+” indicates that the algorithm performed excellently compared to others, while “−” indicates poor performance. If there is no significant statistical difference, it is marked as “=”. The bold text indicates that the algorithm performed the best compared to others. From the data in Table 2, we observe that ICMOCSO performed slightly less well on ZXH_CF4, ZXH_CF7, ZXH_CF8, ZXH_CF10, and ZXH_CF13. However, it demonstrated superior performance on the other 11 test functions. This finding confirms the feasibility of the improvements made to CMOCSO and lays a solid foundation for the subsequent design of multi-PID controller parameter optimization algorithms.

Simulation setup

The parameter range for each individual in the constructed constrained multi-objective problem is 12-dimensional: $K_{p1},K_{p2},K_{p3},K_{p4}\in [22,83.6]$, $K_{i1},K_{i2},K_{i3},K_{i4}\in [0,894]$, $K_{d1},K_{d2},K_{d3},K_{d4}\in [-0.8,0.8]$. Simulation experiment parameters: simulation time is $1s$, population size is 100, maximum iterations are 10,000, maximum overshoot ( ${\sigma }_{max}$) is 1, lower bound of gain margin ( $p_{min}$) is 6, the upper bound of phase margin ( $g_{max}$) is 60, and the lower bound of phase margin $(p_{max}$) is 30.

Result analysis

The ICMOCSO algorithm can obtain the optimal solution set with a single run. From the optimal solution set, we selected three representative solutions for analysis, as shown in Table 3. In Table 3. the first solution set focuses on minimizing synchronization error, the second set on minimizing overshoot, and the third set on minimizing rise time. The bold indicates the best performance in the respective criterion compared to other solutions.

According to Table 3, we can observe that the three sets of solutions exhibit different characteristics in terms of performance. Specifically, the first solution set performs best in terms of synchronization error and has a better steady-state error compared to the other two sets, though other performance metrics are inferior. The second set of solutions has the smallest overshoot, but its other performance metrics are relatively poor. The third solution set has the shortest rise time, and except for steady-state error, its other performance metrics are superior to the other two sets. This finding reveals that in the process of constrained multi-objective optimization of multi-PID controller parameters, we can achieve diverse control effects to meet various control requirements. For different control requirements, we can flexibly select and adjust the parameters of the multi-PID controller to achieve the optimal control effect.

Comparison with other methods

The model constructed in this article is a constrained multi-objective model. Compared with other multi-objective algorithms, constraints are transformed into two objective functions, $f_{less}$ and $f_{greater}$. For any individual $i$ in the population, the constraints are handled using the following method.

(16) $$f_{greater}(i)=\frac{GM(i)-g_{min}}{max(GM(i)-g_{min})}+\frac{PM(i)-p_{min}}{max(PM(i)-p_{min})}$$

(17) $$f_{less}(i)=-\frac{PM(i)-p_{max}}{max(PM(i)-p_{max})}+\frac{{\sigma }^c(i)-{\sigma }_{max}}{max({\sigma }^c(i)-{\sigma }_{max})}$$

In Eqs. (14) and (15), GM represents the gain margin, PM represents the phase margin, $g_{min},p_{max},p_{min}$ are their respective bounds, ${\sigma }^c$ denotes the overshoot, and ${\sigma }_{max}$ is the upper limit of the overshoot.

After processing the constraints, the constructed constrained multi-objective model is transformed into a multi-objective optimization model with nine objectives. Therefore, we chose the upgraded NSGA-III algorithm, suitable for many-objective optimization, as a comparison algorithm.

The method proposed in this article is compared with other PID parameter optimization methods, including the Ziegler-Nichols (Z-N) (Ziegler & Nichols, 1942), the PID parameter tuning method using MATLAB’s built-in Rltool toolbox (Zhang et al., 2016), the NSGA-III method (Deb & Jain, 2013), the GFMMOEA method (Tian et al., 2018), the MaOEAIT method (Sun et al., 2019), and the CMOCSO method (Ming et al., 2022a). In the traditional PID method and the Rltool toolbox method, the master and slave cylinder controllers use the same PID parameters. Among the optimal PID parameter solutions obtained by other methods, to ensure a fair and comprehensive comparison, we selected a solution that performs relatively well in various performance metrics (i.e., a knee point Zhang, Tian & Jin, 2014) for comparative experiments.

The detailed PID controller parameters and their performance are shown in Table 4. In the Table 4, “+” indicates that the optimization method is superior to the other comparison methods in a particular performance metric, while “−” indicates that the method is inferior.

From the table data, we can observe that the optimization method proposed in this article improves performance metrics across the board compared to the traditional PID method, except for overshoot. Compared to the Rltool method, the ICMOCSO method demonstrates better performance in synchronization error, overshoot, settlement time, and integral time square error. Compared to the NSGA-III method, the ICMOCSO method shows advantages in synchronization error, overshoot, steady-state error, and integral time square error. Against the GFMMOEA method, the ICMOCSO method excels in synchronization error, overshoot, rise time, and integral time square error. Compared to the CMOCSO method, the ICMOCSO method outperforms in synchronization error, overshoot, settling time, and integral time square error. Finally, compared to the MaOEAIT method, the ICMOCSO method shows better results in synchronization error, overshoot, settling time, and integral time square error. In summary, the comparative experiments indicate that, on average, the ICMOCSO method reduces the synchronization error by approximately 23.28% and the integral time square error by approximately 46.21% compared to other PID parameter optimization methods.

When constructing the constrained multi-objective problem model for the multi-PID controller, frequency domain indicators were considered to ensure that the optimized PID controller parameters possess anti-interference capability. To test this anti-interference capability, white noise with a noise power of 0.002 and a sampling time of 0.1s was added to the system. The white noise is shown in Fig. 5.

Incorporating the white noise into the simulation model, we again observed the step response curves of each subsystem, as shown in Fig. 6. In Fig. 6, it is clear that under noise interference with a power of 0.002, the multi-PID controller subsystems optimized by the ICMOCSO method, except for the main structure’s PID controller subsystem, reached and maintained within the final steady-state value’s 2% error band the fastest. This demonstrates that the PID parameters optimized by the ICMOCSO method have excellent anti-interference capability, significantly reducing the impact of external disturbances on the system.