Ant colony optimization-based adjusted PID parameters: a proposed method

Preliminary

The ant colony algorithm (ACA) is a new heuristic algorithm proposed (Colorni, Dorigo & Maniezzo, 1991), which can reasonably simulate the foraging behavior of ants in nature. ACA does not need to describe the problem and has the characteristics of global optimization, internal parallelism, positive feedback, and robustness. Besides, compared with other algorithms, the ACA has higher reliability, faster, robust searchability, and easy-to-use implementation in solving combinatorial optimization problems, so it has attracted diverse attention that results in wide use in various fields. Moreover, ACA currently has become a hot topic in the field of multi-disciplines. For example, it is implemented to solve the traveling salesman problem (TSP) with good results. The existing research shows that it can effectively solve some NP-class problems with strong computational power, such as quadratic assignment, vehicle path planning, and job shop scheduling since its enhanced learning system with the characteristics of distributed computing, strong robustness, and easy-to-integrate structure can fit into other optimization algorithms. However, it has some problems, such as the long search time and appearing to be a pause phenomenon.

The implementation of the ant colony algorithm

The ACA is composed of components as follows:

(1) Node establishment: In the upper node, the ant foraging is a process from the starting point to the halfway node and finally to the food source node. Therefore, the practical problems studied should be analogically abstracted with ant colony foraging. This is the first step to establishing the algorithm, and the node establishment itself is the modeling of the actual problem, which is related to the feasibility and advantages of the algorithm. Unreasonable node establishment will lead to the actual problem not being solved and result in low efficiency. Therefore, it is necessary to collect relevant information and establish influential nodes for practical problems in the link of node establishment.

(2) Colony size and total foraging times: “Ant colony” generally expresses the number of ants and describes the total number of foraging times of the entire ant colony. However, due to the limited characteristics of the algorithm, the specific number of ants must be given in practical applications to represent the scale of the ant colony and the total number of foraging times of the entire ant colony. The size of the ant colony, showing the number of ants and the total number of foraging times of the ant colony showing the number of iterations, are extremely important parameters to the algorithm. If the number of ants and iterations is picked too small, the algorithm may be inefficient and eventually lead to the termination of the run before the algorithm convergence. Thus, the result would not be the optimal solution. On the contrary, if the number of ants and iterations is chosen too large, the algorithm will eventually converge. Still, subsequent calculations will be redundant since the number of iterations required for algorithm convergence is exceeded, which will lead to a very long execution time and eventually lose the advantages of ACA. Therefore, when the algorithm is run, it is necessary to set the appropriate number of ants and iterations based on the problem. Additionally, timely adjustments can be made according to the execution effect of the algorithm.

(3) Path evaluation criteria: In the process of ant foraging, the quality of the path can be intuitively expressed by the length of the path, which is an evaluation criterion of the path. However, in the practical application of the ACA, it may be necessary to specify a reasonable index according to the actual demand to replace the path length, which is the formulation used for evaluating path evaluation criteria. The path evaluation criterion is often closely related to the actual intent and represents the definition of the good or bad solution obtained by the final algorithm. Therefore, the formulation of the path evaluation criteria is very related to the final solution of the algorithm representing the solution state of the actual problem.

(4) Transfer probability rule: Transfer probability is an essential basis for ants to transfer from the current node to another node, which is expressed by (2) $$p_{\text{ij}}^k(\text{t})=\{\begin{array}{l}\frac{{\left[{\tau }_{ij}\left(t\right)\right]}^{\alpha }.{\left[{\eta }_{ij}\left(t\right)\right]}^{\beta }}{{\displaystyle \sum _{s\in J(i)}{\left[{\tau }_{ij}\left(t\right)\right]}^{\alpha }.{\left[{\eta }_{ij}\left(t\right)\right]}^{\beta }}},\text{j}\in J_i\\ 0,\text{j}\notin {\text{J}}_{\text{i}}\end{array}$$

${\mathrm{p}}_{\mathrm{i}\mathrm{j}}^{\mathrm{k}}$ represents the probability of ants moving from node i to node j when the number of ant colony iterations equals k generation, ${\tau }_{\mathrm{i}\mathrm{j}}$ denotes the pheromone concentration from node i to node j, ${}^{\eta }{}_{\mathrm{i}\mathrm{j}}$ represents the heuristic function from node i to node j, which can be set according to the actual performance requirements, α designates the weight factor of the pheromone, and β represents the weight factor of the heuristic function. The weight factor represents the degree of influence of each factor on the transition probability. When α > β, the transfer probability mainly depends on the pheromone concentration. When α < β, the transfer probability depends primarily on the heuristic provided that α and β are greater than 1. The transfer probability of the current node to each adjacent node is equal to the ratio of the excitation of each alternative path time, the concentration of the pheromone, and the excitation of the total alternative path times, the attention of the pheromone. Figure 2 depicts it.

Figure 2 depicts that the transfer probability of node 0 to node one is computed by

${\mathrm{p}}_{01}={\displaystyle \frac{\left({1.5}^{\wedge }2\right)\ast \left({1.5}^{\wedge }3\right)}{\left({1.5}^{\wedge }2\right)\ast \left({1.5}^{\wedge }3\right)+\left({1.5}^{\wedge }2\right)\ast \left({1.5}^{\wedge }3\right)+\left({1.5}^{\wedge }2\right)\left({1.5}^{\wedge }3\right)+\left({1.5}^{\wedge }2\right)\ast \left({1.5}^{\wedge }3\right)+\left({1.5}^{\wedge }2\right)\ast \left({1.5}^{\wedge }3\right)}=0.2}$

Similarly, the transfer probabilities of other paths can be calculated and obtained as follows: ${\mathrm{p}}_{02}$ = 0.2, ${\mathrm{p}}_{03}$ = 0.2, ${\mathrm{p}}_{04}$ = 0.2, ${\mathrm{p}}_{05}$ = 0.2. When the transfer probability is obtained, the ant’s movement has a basis, and the algorithm establishes certain rules.

(5) Pheromone update rule: The pheromone update rule directly affects the algorithm’s convergence. If the pheromone concentration increases too fast, the ants will not have time to search for other paths and nodes. Due to the high proportion of pheromone growth in some paths, the ants will have gathered at a part of the paths and nodes, and the final optimal solution obtained by the algorithm would not be a globally optimal solution. Figure 3 depicts it.

Figure 3 depicts that the activity range of the ant colony is limited within the blue coil, which leads to left nodes and paths that are not searched. Although the “optimal solution” can be obtained by premature local search, the optimal herein is not a real optimal in the real sense, but a local optimal that is a pathological manifestation of the ACA because it is not certain that the path formed by the combination of missing nodes and local nodes must be worse than the local optimal solution. Therefore, a reasonable pheromone growth rule must be adopted when designing the updating rule of information number to avoid the algorithm falling into the local optimum prematurely and causing the loss of the solution.

The architecture of the proposed method

The first step is to set the ACA to adjust the PID parameters as the experimental group. MATLAB 7.9 version is used to write the program and three parameters, Kp, Ki, and Kd, are obtained by the path formed by each ant in each generation, and then Kp, Ki, and Kd are loaded into the PID simulation model written in the Simulink in MATLAB 7.9. Afterward, a simulation is run to import the output data of the PID simulation system into MATLAB 7.9 through the output module. Then, the output data is analyzed to obtain the performance value of the PID system formed by Kp, Ki, and Kd values at this time. The next step is to update the pheromone with the set pheromone rules. The steps above are repeated until the algorithm is complete and the relevant data is recorded.

The second step is to set the engineering tuning method to adjust the PID parameters as the control group. The PID simulation model uses the engineering tuning method to calculate the values of Kp, Ki, and Kd, and then the simulation output formed by these three parameters is imported into the program, and the data is analyzed by the same method in the experimental group to obtain the relevant performance values.

The third step is to draw the output waveform diagram of the two groups of experiments, compare and analyze the waveform and related performance values, and draw a conclusion. Figure 4 depicts it.

PID simulation system model

The purpose of building a PID simulation system with Simulink is to simulate the most real PID application effect through an effective model. Second, this article needs to analyze the data through the simulation output to get the pros and cons of the relevant performance index judgment algorithm.

A controlled system function is defined by G(s) = ${\displaystyle \frac{130}{s^2+25s}}$, which can be set according to the actual situation. Figure 5 shows that the system architecture is roughly the same as the PID schematic diagram. This model passes the PID control system through the step input (target temperature) and then acts on the controlled system. Finally, the output of the controlled system is divided into three ways: one is fed back to the PID control system, the other is connected to the oscilloscope (display waveform), and the final is connected to the output module (output data to the work area).

ACA-based PID parameters tuning: experimental group

The construction of the experimental group uses the implementation steps of the ACA. The main steps are given as follows: The first step is the node and path establishment. Kp, Ki, and Kd are found and transformed into a node and path problem. In this experiment, Kp, Ki, and Kd are set as nodes on different axes of different three-dimensional coordinate axes.

Kp, Ki, and Kd correspond to the x-axis, the y-axis, and the z-axis, respectively, and are set with 50 nodes from 0.2 to 10 with a step value of 0.2. Figure 6 depicts that a line of the same color and a node connection represents a valid path, and the sequence is the Kp connected to Ki, and Ki is then related to Kd. Therefore, the rules of nodes and paths are established, and when the ant selects a node on a different axis, the ant’s Kp, Ki, and Kd values are obtained.

The second step determines ant colony size and total foraging times. In the program, the number of ants (count) represents the colony size, the number of iterations (N.C.) means the total foraging times and the count is set to 50. N.C. is set to 200; the number of ants in one generation is 50, and the total number is 200 generations. At the end of every 50 ants searching the path, the iteration counter is increased by one until the number of iterations reaches 200 and the entire program is finalized.

The third step is path evaluation criteria. In this experiment, the comprehensive performance value (CPV) is set to replace the path length as the path evaluation criterion. First, the performance indicators related to step response are introduced (the output of the controlled object of the simulation system is step response).

• (1) Overshoot (MP): the maximum value of the simulation result.

• (2) Rise time (tr): the time from the start to the first time for the specified stable value.

• (3) Adjustment time (ts): the time of transient process experienced when the deviation between the actual output and the specified stable value reaches the allowable range (1% in this experiment) (timing from t = 0)

• (4) Steady-state error (err): specifies the difference between the stable value and the steady-state output state.

Comprehensive performance value (CPV) calculation formula: (3) $$\mathrm{C}\mathrm{P}\mathrm{V}={\displaystyle \frac{{10}^{32}}{M{P}^{20}\ast t{r}^3\ast t{s}^3\ast er{r}^4}}$$

Equation (3) implies that the higher the CPV, the better the performance of the control system, and the CPV is more biased to the performance value of overshoot (M.P.).

Step 4: Transfer probability rules and transfer implementation.

We mainly establish the correlation function and parameters. To pursue the speed of the control system, the enlightening function in this design is set as follows:

• (1) Path heurism function from start point to Kp node: ŋ = 1 + 0.1 × Kp (min = 1.02, max = 2)

• (2) Path heurism function from Kp node to Ki node: ŋ = 1 + 0.1 × Ki (min = 1.02, max = 2)

• (3) Path heuristic function from Ki node to Kd node: ŋ = 1 + 0.2/Kd (min = 1.02, max = 2)

The pheromone concentration function is set as follows: the pheromone concentration of each path is initialized to 1 and changes as the pheromone concentration is updated. The greater the inspiration is equivalent to the larger the Kp value, the larger the Ki value, the smaller the Kd value, the faster the system recovery, which can effectively reduce the rise time and adjustment time, but the corresponding overshoot and stability error will increase, which is equivalent to a trend of one part of the performance, away from the other part of the performance. This setting will make the control system faster. However, stability is not guaranteed, so the weight of the inspiration function should not be too large and less than the weight of the pheromone concentration function because the inspiration is very limited. The pheromone is focused on the global. The hope is to focus on the control system’s speed to ensure comprehensive performance rather than mindlessly pursuing the control system’s speed. Therefore, in this design, the pheromone concentration weight α is set to 10, and the heurizing weight β is set to 1.

The fifth step is the pheromone update rule: the update frequency is set to tune the algorithm once every iteration occurs, and each generation’s optimal path (the path with the largest CPV) is selected for update. To prevent local search, the update of pheromone concentration should not be too large, so it is set to increase CPV/100 each time.

At this point, the experimental group was completed. Since K.P., K.I., and K.D. are subdivided into 50 pieces of data from 0.2 to 10 with a step value of 0.2, the composed path reaches 50 × 50 × 50 = 125,000 kinds, and the amount of data is too large to be displayed. The algorithm is a dynamic process that cannot exhaust each state’s specific state. Therefore, the following K.P., K.I., and K.D. are intercepted from 0.2 to 1, and the initial stage of the algorithm is used as a program flow demonstration.

Adjusted PID parameters by engineering setting method: control group

The methods of tuning PID parameters in engineering mainly include critical proportion, response curve, and attenuation. This article uses the 4:1 attenuation curve method to adjust PID parameters to construct a control group. The steps are given as follows:

(1) The regulator integral time is set to infinity, the differential time is set to zero (Ki = ∞, Kd = 0), Kp is appropriate, and the control system is input by pure Kp action. After stability, the proper proportion is reduced until the adjustment process changes reach the specified 4:1 attenuation ratio. The relevant parameters of proportionality δs and attenuation operation period T.S. in the case of 4:1 attenuation are obtained.

(2) According to the equation, critical scale δs and the calculation of the critical period Tk is defined by

Kp = 1/delta s; Ki = Kp/(0.3 * TS); Kd = Kp * (0.1 * TS).

(3) According to the actual effect correction.

It is measured that when δs = 0.37 (Kp value in step 1), the two adjacent peak values are 74.73 and 56.22, respectively, and the peak time is 1.057 and 1.167 s, respectively. (74.73 − 50)/(56.22 − 50) = 3.96, which is close to the 4:1 attenuation ratio, so it meets the conditions. Then T.S. = 1.207 − 1.057 = 0.15.

Kp = 2.7, Ki = 0.045 and Kd = 0.015 can be obtained from the step equation and final correction.