Comparison of predictive control using Self-Organizing Migrating Algorithm and MATLAB fmincon function

The aim of this paper is to evaluate the usability of the self-organizing migrating algorithm (SOMA) in a nonlinear system predictive control area. The model predictive control is based on an objective function minimization. Two approaches to model predictive control applied on a nonlinear system are studied here. Firstly, the SOMA was used to minimize the objective function, secondly, the fmicon function included in the MATLAB optimization toolbox was used for the same. The nonlinear system simulated here is an exothermic semi-batch reactor mathematical model based on a real chemical exothermic process. Also the input data used here to simulate the process were obtained from the same real process. Results obtained by the simulation means were than evaluated using suitable criterion which was defined for that purpose and discussed.


Introduction
Control of nonlinear systems brings challenges in the controller design.The current availability of powerful computing technologies enables using of complex computational methods.One of such complex method is also the self-organizing migrating algorithm (SOMA).This algorithm can be used for solving of various optimization problems.Such problem definitely is also the model predictive control (MPC).Here suitable algorithm minimizes an objective function which is based on the responses from a real system model and the real system itself.Minimizing the objective function using SOMA is studied here and the comparing with MATLAB fmincon function is also done to evaluate the SOMA control ability.
The model of real process on which are the simulations performed comes from leather waste processing procedure.Part of the recycling process includes also an oxidation-reduction reaction which is strongly exothermic and can be controlled by the chromium sludge (the waste produced by leather industry) into the hot reaction blend of chromium sulphate acid dosing [1].

The nonlinear system to be controlled 2.1 Semi-batch process
As was already mentioned, the nonlinear system here represents the exothermic semi-batch reactor in which the chromium leather waste is recycled.The chemical reactor is a vessel with a double wall filed with a cooling medium.It has a filling opening, a discharge outlet, cooling medium openings and a stirrer.
The reactor is filled with initial filling given by the solution of chemicals without the chromium sludge (filter cake).The sludge is fed into the reactor to control the developing heat since the temperature has to stay under a certain critical level (T(t) < 373.15K), otherwise the reactor could be destroyed.On the other hand, it is desirable to utilize the maximum capacity of the reactor to process the maximum amount of waste in the shortest possible time (higher temperature is desirable).Therefore, an optimal control strategy has to find a tradeoff between these opposite requirements.

System mathematical model
Based on the balanced equations (the mass and heat balance), system mathematical model was derived [2].The equations describing the system are displayed here (Eq.1-4): Individual symbols have the following meaning: m is the total weight of reaction components in the reactor, a is the mass concentration of the reaction component in the reactor, c = 4500 J•kg•K -1 is the specific heat capacity of the reactor content and T its temperature.F I , T I = 293.15K and c I = 4400 J•kg•K -1 are the reaction component input mass flow rate, temperature and specific heat capacity. 1 and m C = 220 kg are the cooling water mass flow rate, input temperature, output temperature, specific heat capacity and weight of the cooling water in the cooling system of the reactor, respectively.
Other constants: The fed-batch reactor use jacket cooling, but the effective heat-transfer area (S = 7.36 m 2 ) in the mathematical model was treated as constant, not time varying.The initial amount of material placed in the reactor takes about two-thirds of the in-reactor volume and the reactor is treated as ideally stirred, so we can do this simplification.
Variables F I , F C , T I , T CI , can serve as manipulated signals.However, from practical point of view, only F I and F C are usable.The T I or T CI temperature change is inconvenient due to the economic reasons (great energy demands).

System constrains and other limits
Maximum filling of the reactor is limited by its volume to m = 2450 kg approximately.The process of the chromium sludge feeding F I has to be stopped by this value.Practically, the feeding F I can vary in the range 3 ; 0 ∈ I F kg.s -1 .As stated in the system description, the temperature T(t) must not exceed the limit 373.15 K; this tem-perature value holds also for the coolant (water) but it is not so critical in this case as shown by further experiments.

Control of the system
Our nonlinear system here represents a chemical reactor.The state of art of chemical reactors control presents for example Luyben in [3] and [4], control and monitoring of batch reactors also describes Caccavale et al. in [5].Generally, it can be stated that chemical reactors controllers use various control methods, such as PI controllers, adaptive control methods, robust approaches, predictive control and the like [6][7][8][9][10][11][12][13][14][15].The model predictive control [16][17][18] belongs to the one of the most popular and successful approaches for semi-batch reactors control.However, this methodology brings some difficulties in finding optimal control sequence especially when complex nonlinear model is utilized.Interesting way how to cope with the optimization problem offers the usage of evolutionary algorithms [19][20].Some review of the recent state can be found in for example in [21]

Model predictive control
Two different approaches to the model predictive control of the given system are introduced in this paper.At first, the model predictive controller uses SOMA algorithm for the optimization of the control sequence.This methodology ensues from model predictive control method [22] while it uses same value of the control signal for whole control horizon in order to reduce computational demands of the controller.Secondly, the classic MPC controller, which uses Matlab Optimization Toolbox fmincon function, was used.
The main idea of MPC algorithms is to use a dynamical model of process to predict the effect of future control actions on the output of the process.Hence, the controller calculates the control input that will optimize the performance criterion J (Eq. 5) over a specified future time horizon [23]: where k is discrete time step, N 1 , N 2 and N u define horizons over which the tracking error and the control increments are evaluated.The u t variable is the tentative control signal, y r is the desired response and y ˆis the network model response.The parameters λ and ρ determine the contribution that the sums of the squares of the future control errors and control increments have on the performance index.Typically, the receding horizon principle is implemented, which means that after the computation of optimal control sequence only the first control action is implemented.Then, the horizon is shifted forward one sampling instant and the optimization is again restarted with new information from measurements.

SOMA algorithm
The Self-Organizing Migrating Algorithm (SOMA) is based on the self-organizing behavior of groups of individuals in a "social environment".It can be classified in two ways -as an evolutionary algorithm or as a so called memetic algorithm.
SOMA algorithm can be used for optimizing any problem which can be described by an objective function.This algorithm optimizes a problem by iteratively trying to improve a candidate solution, i.e. a MATEC Web of Conferences 210, 02041 (2018) https://doi.org/10.1051/matecconf/201821002041CSCC 2018 possible solution to the given problem.The SOMA has been successfully utilized in many applications [24][25][26], while interesting comparison to with simulated annealing and differential evolution is provided by Nolle et al. in [27].

SOMA simulations
Simulations were performed in the Mathematica 8.0 software.Here the algorithm SOMA was used for the cost function (5) minimization and was set as follows: Migrations = 25; AcceptedError = 0.1; NP = 20; Mass = 3; Step = 0.3; PRT = 0.1; Specimen = {0.0,3.0, 0.0}; Algorithm strategy was chosen All To One.First two parameters serve for the algorithm ending.Parameter "Migrations" determines the number of migration loops, "AcceptedError" is the difference between the best and the worst individuals (algorithm accuracy).If the loops exceed the number set in "Migrations" or "AcceptedError" is larger than the difference between the best and the worst individuals, the algorithm stops.Other parameters influence the quality of the algorithm running."NP" is the number of individuals in the population (its higher value implicates higher demands on computer hardware and can be set by user), "Mass" is the individual distance from the start point, "Step" is the step which uses the individual during the algorithm, "PRT" is a perturbation which is similar to hybridizing constant known from genetic algorithms or differential evolutions."Specimen" is the definition of an exemplary individual for whole population.For details see [28].Seven different simulations using SOMA algorithm were performed.First three simulations (SOMA1 -SOMA3) were done to study the control horizon N u influence, next three (SOMA4 -SOMA6) the prediction horizon N 2 influence and the last one (SOMA7) is the simulation with an optimal setting.All settings can be seen in Table 1.The control horizon (N u ) actually means the time interval, for which the actuating variable (F I ) has constant value.It is generally better to set it as short as possible because of more rapid influence on the system, but on the other hand the lower value increases the computing time during the calculations.So it is necessary to find the control horizon value, which balance between these two requirements.
The prediction horizon (N 2 ) determines how forward controller knows the system behavior.If the horizon is too short, the controller doesn't react in time and the system may become uncontrollable.Long horizon means again the more demanding computation, i.e. the need of more powerful computer hardware.
Graphical output of SOMA7 (the optimal settings) simulation is depicted in Figure 1.The two most important dependencies are here -the in-reactor temperature and the chromium sludge dosing development.As was already mentioned, the temperature has to stay under critical point 373.15K.The chromium sludge dosing shouldn't embody any rapid changes.

Conventional MPC approach
This part was simulated using Matlab/Simulink, where the standard Matlab Optimization Toolbox function fmincon with receding control strategy was implemented.The fmincon function used trust-regionreflective algorithm [29].
To get the similar settings as in the SOMA case (the constant control action for the whole length of the control horizon Nu= 60), the sample time was set to 60s.The control horizon N u and the prediction horizon N 2 were set to 10.The rest of the controller design remained same -the predictor was based on the white-box model described by equations (1 -4), cost function used by MATLAB was also the same (5).
The first set of simulations showed that the control did not provide acceptable results.The permanent control error and/or controlled variable overshoot where not suitable here.It was found that the problems were located mainly in the beginning of the control process.The controller took enormous control actions there.This strange behaviour was result of reaction kinetics and strongly exothermic reaction combination.Even small concentration growth of the chromium sludge (the increase in actuating variable) causes steep rise of the temperature, but the reaction kinetics can cause a response delay to the dosing.
To prevent this unwanted behaviour, new criterion based on the criterion (5) was defined.This enhanced criterion was able to penalize values of the control signal in the process start part.Also, at the same time the penalization has to decrease taperingly.The new https://doi.org/10.1051/matecconf/201821002041CSCC 2018 enhanced criterion is described by equations (6-7) and the controller settings are placed in table 2.
The parameter γ c defines here the speed of the decrement in γ.In this way we can influence the speed of dosing the chromium sludge, the actuating value.We can say that γ parameter defines penalization of the control signal, while the ratio γ/γ c specifies the length of the penalization interval.Too high γ parameter or γ/γ c ratio caused delays or oscillations (the settings MLB2 in Table 2).On the other hand, small γ/γ c ratio led to overshoots of the temperature (the settings MLB3 in Table 2).The best result obtained using this approach was obtained for MLB4 settings and is displayed in Figure 2.

Results comparison
Results of the best SOMA and MATLAB simulations were selected for the comparison.To compare the control error, the criterion function S y was defined (Eq.8): Other criterion S u (Eq.9) was defined to monitor the speed of the control signal changes.From the practical view, the monitoring of it is very important, because lifetime of the mud pump (actuator) that injects the chromium sludge to the reactor would be shortened significantly in case of steep changes.
The number of steps computed for the criterions S y and S u is defined by t f and was set to 50 steps.
Observed were also the maximum overshoot of the output value y max and the time of the reaction (dosing) t b .
For the reason that the plant is very exothermic and it is very sensitive to the exceeding of the desired value of the temperature (y r = 370K), it was necessary to observe the maximum overshoot of the output value y max .Furthermore, it is essential to observe the time of the reaction (chromium sludge dosing) t b .The heating up and maintaining the system temperature usually takes about 3000s and after that only cooling is performed.
In fact, there was not significant difference in the temperature overshoots between the SOMA and MATLAB, they were quite similar.Anyway, as can be seen in Table 3 the result obtained by SOMA was a bit better.Also the results provided by criterion S y were in both cases close.The lower value is better value in case of S y and again the SOMA control quality prevailed.The time of dosing achieved by SOMA was shorter approximately for one minute (58 seconds).
On the other hand, the MATLAB gave better results for the S u criterion.The SOMA value 2.3200 was higher than the 1.5500 MATLAB value.The actuating device would last longer without servicing in MATLAB case.

Conclusion
Although the SOMA and MATLAB were both able to control the complex nonlinear process, there were some differences.A surprising difference emerged when the MATLAB was not able to provide satisfactory control results using the same objective function as SOMA algorithm did.That is why the purpose function had to be changed for the MATLAB simulations.After this change, the comparison was made from four points of view.Firstly, the value of the in-reactor temperature overshoot and the related quality of the in-MATEC Web of Conferences 210, 02041 (2018) https://doi.org/10.1051/matecconf/201821002041CSCC 2018 reactor temperature course were observed.Secondly, the time of processing which is important for effectiveness of a real plant and also the course of the actuating signal that is important from the practical point of view were monitored.Although the results were similar, SOMA showed generally better results than MATLAB.
The results show that the nonlinear system control can be successfully improved by evolutionary algorithms.