Control of Linear Multichannel Objects with Numerical Optimization

1. Introduction

Control of objects with two inputs and two outputs, where each input affects each output, requires careful design of the controller [1,2,3,4,5,6,7,8,9,10]. In particular, examples of such objects are the problem of simultaneous stabilization of the power and frequency of a semiconductor laser by controlling its current and temperature, with the temperature controlled by a Peltier- based microcooler. Another example of such problems is the automatic mixing of several substances, when it is necessary to simultaneously control both the increment of the mixture volume and the concentration of each substance. Another example of such a problem is the well-known problem of filling a swimming pool, with the need for separate control of the flow rate and temperature of the water. Similar problems arise in controlling the movement of any vehicle, especially space, air and water. There are many examples of such problems in robotics, chemistry, metallurgy, and scientific research.

Let us consider the problem of controlling linear objects. The proposed approach can also be applied to nonlinear objects, but the obtained result should be more carefully checked for the case of test signals of different amplitudes. In general, the problem of controlling nonlinear objects is much more complex; it is not possible to consider it thoroughly enough within the framework of this article, except for such types of nonlinearity as limitation or relatively small deviation of the characteristic from the linear one (no more than 10%). The discussed method is unsuitable for controlling objects with nonlinearity of the “dead zone” type, as well as “dry friction”, “hysteresis”, or other ambiguous nonlinearities. The proposed method is insufficient for controlling non-stationary objects; adaptive control systems are required for such objects.

The mathematical model of a linear two-channel object is defined in the form of a matrix transfer function [11,12].

$$W\left(s\right)=\left[\begin{array}{cc}{w}_{11}\left(s\right)& w_{12}\left(s\right)\\ w_{21}\left(s\right)& w_{22}\left(s\right)\end{array}\right].$$

(1)

There are many articles that talk about an object of size n×n, but no articles have been found with reliable information about successful control in such a formulation of the problem, i.e. with multichannel negative feedback and a multichannel sequential controller, even for objects of size only 3x3. Therefore, the term “multichannel objects” should be recommended to be used more carefully than is currently accepted in the literature. “Two-channel” would be more correct if even “three-channel” ones are not involved. However, we have verified our theses for validity by induction on the example of an object of size 3x3 and, if necessary, are ready to go further, to objects of size 4 x 4, but we do not know of such problems from practice, so we are in no hurry with these studies.

This approach should not be confused with multiparameter optimization problems or multicriteria optimization problems.

Control in a negative feedback loop is carried out according to the structure shown in Figure 1.

This structure best satisfies the requirement that the output signal should match its assigned value, called “Task”. Here, the element marked with a circle with a sign $“\mathsf{\Sigma }\mathrm{”}$ inside is an adder, but if there is a minus sign before the input, this input is subtractive. This and the following illustrations are made in the VisSim 5.0e software, which has some features for displaying various mathematical symbols, numerical values, and other elements. In particular, the minus sign looks like a hyphen, blue rectangles indicate composite blocks that have a complex structure inside, and their name is set by the user. Arrows indicate the direction of signals. When using the structure according to Figure 1, the task of controlling an object is reduced to calculating a mathematical model of the regulator (controller) based on a known mathematical model of the object.

2. Statement of the Problem

In the early stages of the development of automatic control theory, the control problem was posed and solved in terms of finding a control signal that must be fed to the input of the object so that its output values (called output signals) coincide with the prescription “Task”. This approach is impractical. The task can change, therefore, calculating the control signal loses its relevance. In addition, the object is always affected by factors that cannot be measured and taken into account in such a problem statement. The use of the structure according to Figure 1 is a different problem statement: instead of calculating the signal controlling the object, a mathematical model of the device is calculated, which itself calculates the necessary control signal depending on the task and the actual output signal of the object. A characteristic property of systems with negative feedback is that with a correct calculation of the regulator, the control problem is solved successfully even in the presence of unknown and uncontrolled effects on the object, which change its output signal in an unpredictable way. In the frequency band where the transfer function of the controller and the object connected in series with it is much greater than one, the output signal coincides with the task quite accurately. In this case, the difference between the output signal and the task, called the error, is inversely proportional to the product of the controller transfer function and the object transfer function. The problem with using this method is that any object has a transfer function limited by the frequency band. With increasing frequency, the transfer function decreases, although in some local frequency ranges it can increase. However, starting from some frequencies, it decreases so sharply that no amplification of the input control signal can achieve a significant response at the object output. For this reason, the transfer function of the entire loop in a conditionally open form inevitably reaches a value close to unity, and then decreases even more. In the frequency range where this value is much less than unit, the control loop no longer has a noticeable effect on the object, so its output signal in this frequency range is not controllable. We have to put up with this. But in the intermediate frequency band between those where the transfer function of the circuit is much greater than unit, which allows the object to be controlled, and those where the transfer function is much less than unit, which excludes the possibility of controlling the object, there is a frequency region where the transfer function of the circuit lies in the range from 0.1 to 10, i.e. from $-20dB/dec$to $+20dB/dec$. In this frequency range, the behavior of the transfer function of the loop affects the stability of the system. If at these frequencies the signal delay exceeds $180°$ (or $\pi rad$ which is the same value in other measure units) then the system will be unstable. Therefore, the art of the designer of the locked dynamic systems consists in calculating such a controller that would ensure the stability of the locked system, as well as its sufficient accuracy in the frequency region where control is provided.

For multichannel systems, this is not enough. An additional property called “autonomization” or diagonal decoupling is required. It consists in the fact that any change in the task signal at any input should not only provide a corresponding change in the output signal at the corresponding output, but also that it should not affect the output signals of all other inputs. This requirement can be met only approximately. There was a scientific direction called “invariant control” that studied such control methods so that the control error was always strictly zero, but this area of automatic control theory has lost its relevance. The object is always affected by uncontrolled disturbing factors, which are usually described as an unknown interference, which is applied to the output of the object in the form of an additive unknown signal. Instead of the unrealistic task of completely eliminating the influence of disturbance, the modern section of control theory dealing with high-precision systems considers the problem of limiting the error to a certain value measured as a percentage of the task increment. An error of less than 5% of the task was considered acceptable at the early stage of automatic control theory; in modern high-precision systems, requirements are imposed to suppress the error to values of no more than 0.0001% and, in some cases, to an even smaller value. In problems of controlling multichannel objects, as a rule, the error is reduced to values of the order of 1–5% . In this case, we are talking about a dynamic error, since a static error, i.e. an error that is established after a sufficiently long time, provided that the task signals are step functions, can be reduced strictly to zero. This error does not include errors in measuring the output signal, which can never be reduced strictly to zero.

Thus, the problem of controlling an object, taking into account the above, is reduced to the problem of designing a regulator for an object with a known mathematical model using the structure according to Figure 1. Other structures and other formulations of the control problem are also known, which we do not consider in this article.

The controller for a linear object (1) can be represented as a matrix transfer function [13,14,15,16]. For example, for an object with two inputs and two outputs, such a transfer function has the following form:

$$R\left(s\right)=\left[\begin{array}{cc}{r}_{11}\left(s\right)& r_{12}\left(s\right)\\ r_{21}\left(s\right)& r_{22}\left(s\right)\end{array}\right].$$

(2)

In the simplest case, this matrix is diagonal, that is, only the diagonal elements are nonzero, and the remaining elements are equal to zero:

$$R\left(s\right)=\left[\begin{array}{cc}{r}_{11}\left(s\right)& 0\\ 0& r_{22}\left(s\right)\end{array}\right].$$

(3)

This means that there are only feedback connections from each output to the corresponding input and no other connections.

In more complex structures, additional loops may appear, for example, a loop that covers an object in addition to the device that subtracts the output signal from the task. There may also be a loop that covers the controller with feedback. The feedback may contain a transfer function that differs from unity. Ultimately, such modifications can be reduced by equivalent transformations to the structure shown in Figure 1, with the difference that the controller’s transfer function will be different, more complex.

The choice of a controller containing parallel proportional, integrating and derivative channels, which for this reason is called a PID controller, is determined by the fact that for the overwhelming majority of practical problems this choice is sufficient, and also by the fact that simpler structures in which one or two of the specified paths are missing are a special case of such a controller. These simpler structures can be obtained in special result in the same way as the design of a PID controller, which does not change the approach to the problem of controlling an object of type (1).

The development of computer technology and software has simplified the task of designing regulators for complex systems as much as possible [17,18]. For this reason, there is no need to even think about whether the problem can be solved using only a diagonal matrix regulator, or whether a full-fledged multichannel regulator with all non-zero elements of the matrix transfer function is necessary [19,20,21,22,23,24,25,26,27,28]. In many cases, it is sufficient to simply calculate a full PID controller, which in practice for objects of dimension $2\times 2$ and $3\times 3$ turns out to be insignificantly more difficult than calculating simplified diagonal controllers, as will be shown below.

The most common examples in the literature are objects with dimensions $2\times 2$.

In this case, the object model is a matrix transfer function of the same dimension.

For control, a matrix PID controller of the same size is required [13,14,15,16,17,18,19,20,21,22,23,24,25,26], or perhaps a more complex version, for example, a controller containing additional channels of double derivation, double integration, both of these channels, or something more complex structure decision. For a PID controller, its individual elements are described by the sum of the proportional, integrating, and derivative channels with the corresponding coefficients:

$$R_{ij}=k_{Pij}+{\displaystyle \frac{k_{Iij}}{s}}+k_{Dij}s={\displaystyle \frac{k_{Pij}s+k_{Iij}+k_{Dij}{s}^2}{s}}.$$

(4)

The structure of the control system is shown in Figure 1.

3. Theory of Solving Problems of Multichannel Objects Control

The theory of locked dynamic systems is well developed: negative feedback ensures accuracy provided that the entire system is stable. Structural design consists of choosing the structure of the controller; parametric design consists of determining the parameters of the controller based on the chosen structure.

Structural design is not the subject of this article, since PID controllers are the most common, and the reason for this is their simplicity and efficiency. Simpler controllers, such as PI, PD, are not considered, since such simplification does not greatly simplify the method, but limits its capabilities. Complication by introducing additional channels of double integration or (and) double derivation is also not considered, since the approach under consideration can easily be extended to these structures. Structures with a large number of additional loops, such as the Smith predictor, local and pseudo-local connections, etc., are not considered due to their excessive complexity and low efficiency, which has also already been confirmed by modeling and practice.

Parametric design in our case is the main content of the system design task. Methods for calculating the parameters of the regulator can be grouped into large sections: analytical, engineering, tabular, empirical and numerical.

Analytical methods are based on the direct solution of matrix equations, which leads to the Riccati equation. If the matrix transfer function is not degenerate, it can be inverted. Inversion of this matrix is one of the important elements of solving the problem under consideration analytically. Even if the control problem is formulated in a different mathematical formulation, analytical methods cannot do without inverting the matrix transfer function of the object.

An inverse matrix is a matrix that, when multiplied from the right or left by the original matrix, yields the identity matrix. If the matrix transfer function has at least one delay link, the inversion of such a matrix transfer function yields a transfer function that cannot be implemented as any real device.

Let’s look at an example.

$$W\left(s\right)=\left[\begin{array}{cc}{\mathrm{e}}^{-\tau s}& 0\\ 0& {\mathrm{e}}^{-2\tau s}\end{array}\right].$$

(4)

Inverting this matrix gives the following result:

$$W^{-1}\left(s\right)=D\left(s\right)=\left[\begin{array}{cc}{\mathrm{e}}^{\tau s}& 0\\ 0& {\mathrm{e}}^{2\tau s}\end{array}\right].$$

(5)

Indeed, it is easy to verify that multiplication from the left or right by this matrix yields a unit matrix. But the resulting matrix $D\left(s\right)$ consists of functions that cannot be implemented in any physical device. The fact is that the function ${\mathrm{e}}^{-\tau s}$ describes in operator form (in the form of Laplace transforms) such a transformation of the input signal, which consists in delaying the input signal by time $\tau$. However, at the same time, the transfer function ${\mathrm{e}}^{\tau s}$ describes such a transformation, which forms an input signal at the output of the element with a negative delay, that is, $\tau$ seconds earlier than this signal arrived at the input of the element. Such a “signal prediction link” cannot physically exist, it is something like a hypothetical “time machine”. Thus, if there is a delay in the model of the object, then analytical methods cannot be applied to the calculation of the regulator for such an object, or it is necessary to invent methods for reducing the delay to some other mathematical description, for example, to the Padé approximation. But such an approximation makes all calculations excessively complex and yet unreliable, and the result of the calculation is unreliable. The main objection to this approach is based on the fact that any mathematical model of any real object that does not take into account the delays of this object is insufficiently detailed for its effective use. This thesis is substantiated in Appendix A.

Tabular methods such as the Cohen–Kuhn, Ziegler–Nichols, Cohen–Coon, as well as Chien–Hrones–Reswick method and similar ones [29] are always unreliable, so they are hopelessly outdated. These methods are based on the assumption that the mathematical model of the object is described accurately enough by a series connection of a first- or second-order filter and a delay element:

$${W\left(s\right)={\displaystyle \frac{k}{1+2\xi Ts+T^2{s}^2}}\mathrm{e}}^{-\tau s}.$$

(6)

Tabular methods suggest obtaining a graph of the transient process as a reaction to a single step jump, then using this graph to determine several basic parameters of this model. Then, using the table, one should determine the parameters of the PID or PI controller for such an object. A modification of this method is another method, in which the object is covered by proportional negative feedback, after which the coefficient of this feedback is selected so that undamped oscillations of constant amplitude are established in the system. Then, using the frequency and amplitude of these oscillations, the coefficients of the PID controller are also calculated from the table. All such methods rarely give an acceptable result, and they never give the best possible result. The problem with them is that they often result in an unstable system.

Some authors report that after using one of the tabular methods they slightly adjusted the coefficients, as a result of which they obtained satisfactory results. But in this case it should be recognized that the tabular method could not have been used at all; it would have been enough to set arbitrary coefficients of the controller, and then adjust them. This empirical method of adjusting controllers is widely used. It is used much more often than it is stated in articles, because the fact of using the empirical method does not provide grounds for writing a scientific article. The coefficients of the controller are simply successively increased until the stability in the system is violated, after which they are reduced to a value at which stability is restored.

Engineering methods consist of constructing an asymptotic logarithmic amplitude-frequency characteristic, as well as a phase characteristic with subsequent use of the logarithmic Nyquist-Mikhailov stability criterion. These methods work well for single-channel systems, but they have not been developed for multichannel systems, and it is too late to develop them, since they are morally obsolete. Before the active use of computing technology, constructing logarithmic graphs was much more effective and simpler than constructing graphs on a linear scale. An effective technique for constructing such graphs and their use for analyzing the stability of a system was developed, as well as for ensuring stability by correcting the logarithmic amplitude-frequency characteristic. A method for calculating a mathematical model of a correcting device based on the obtained logarithmic amplitude-frequency characteristic of this device is also known. At present, constructing logarithmic graphs that approximately describe the control object is already a more complex technique than numerical optimization. Thus, for linear objects with delay, the best method is a numerical optimization method, and if the object has at least the simplest nonlinear component, then the numerical optimization method is the only reliable method that allows obtaining adequate results with the least labor costs.

The most reasonable approach is a preliminary attempt to design a PID controller [30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58], since its simplification is not advisable, since this only slightly simplifies the design procedure, and its complication may be inappropriate due to the excessive complexity of this procedure and due to the insufficient efficiency of this measure in many practical cases.

The assertion about the insufficient efficiency of more complex structures is based on the fact that many types of complication of regulators lead to non-rough systems. Non-rough systems are those systems in which an insignificant deviation of the values at least one of all coefficients of the mathematical model of the object or regulator causes significant and even fatal changes in the type of transient process, including a violation of the stability of the entire system. Such methods leading to these undesirable results include the Smith predictor method, the method with high-order derivation, the method with a compensator, i.e. an element that is connected in series with the regulator so that in this element the numerator of the transfer function contains completely or partially the problematic factors of the denominator of the transfer function of the controlled object.

The advantage of numerical optimization is the ability to obtain results automatically. Only simple computer is necessary as well as free or inexpensive software and mathematical model of the object, as well as knowledge that can be acquired in the relevant literature, taking into account the recommendations of this article. The disadvantage of the numerical optimization method is the impossibility of solving the problem in general. Any solution is relevant only for a mathematical model with specific numerical values of all parameters of this model. However, some deviations of the actual values of the model parameters from the values used for the calculation are still allowed. If the solution remains relevant with a deviation of these parameters by 1-2%, the solution can be called rough. If the solution remained relevant even with a deviation of the parameters by 5-10%, this solution would be clearly robust. However, this is not always possible. If the model parameters can deviate from the calculated values, then an adaptive control system is required, which is not considered in this article.

Before describing the solution to the problem, we will establish the requirements for the transfer function of the object (1), necessary for the problem to have a solution. First, the matrix transfer function of the object (1) should not be singular. Its rank should coincide with its dimension; in this case, the rank should be equal to 2. In addition, it is highly desirable, although not necessary, that with $s=0$ the matrix $W\left(0\right)$ should also be non-singular. If this matrix is singular, then control in the static mode will have some features [58]. They consist in the fact that even in the static mode, the control signals cannot be static. That is, to ensure that the output signals of the object take the prescribed value in the steady-state mode, it is not enough for the control signals arriving at the input of the object to have some necessary static value; these signals will have to change over time even to ensure a steady-state output state. This can be proven theoretically and experimentally [58].

4. Method for Solving the Problem

In general terms, the solution method consists of using the following initial means: a) software for modeling and optimization; b) test signal; c) modeling conditions (time sampling step and modeling time); d) objective function; e) optimization method; f) result analysis method; and g) method for improving the result if it is not good enough.

The question naturally arises: “If the result is optimal, then how can we raise the question of its further improvement?”

However, the point is that the objective function is not the only option that follows from the control goals. For a set of control goals, the most diverse objective functions can be substantiated, formed and used. In addition, the result also depends on the modeling conditions and on the weight coefficients in these objective functions, on which the method for further improving the result is based.

The optimization method deserves special comments. If the optimization is successfully completed, then it makes no difference which optimization method was used. Any optimization method that gives a result should give the same result within the permissible error. If different optimization methods give different results, this indicates an error in solving the problem. Most likely, this can happen when the objective function has several local extrema. In this case, a designer can distinguish a local extremum from a global one by the value of the objective function at these points. If the minimum of the function is sought, then at the point of the global minimum the objective function takes the smallest value, and at the points of the local minimum, its value is greater.

If one of the numerical optimization methods does not lead to the desired result, and another method does, then the system designer can simply use the results obtained by a more effective method, and there is no need to think about the reasons why the other method was unsuccessful. Therefore, if somebody were able to buy the required product in one store, he is no interested in the reasons for the absence of this product in other stores. In general, various optimization methods can be theoretically analyzed and the different efficiency of various methods can also be theoretically substantiated. Although experimental substantiation is quite sufficient to move on.

Finally, if none of the numerical optimization methods embedded in the software leads to the desired result, it is necessary to find out the reasons for this. Instability of the optimization procedure is most often a consequence of an insufficiently effective objective function. The objective functions and their modifications proposed in this article do not cause problems except for those cases when the problem cannot be solved by numerical optimization methods due to the absence of a global minimum for the mathematical model of the used object. This issue is also discussed below in Appendix A.

In the VisSim 5.0e, as in other versions of this software, there are two types of optimization. One of these types is the search for parameters that give the objective function a zero value. We do not use this type of optimization for the problem under discussion. The other type of the optimization that we use is the search for parameters that give the minimum of the positive-definite objective function, which in this case is called the “cost function” or simply “cost”. We will use this name further.

The method we propose is developed and described in some detail in our publications [11,12,17,18]. To solve the problems of control of multichannel objects, some modification of the cost function is required.

The cost function for this case is the sum of the cost sub-functions for each individual channel [58]. For each channel, such a cost sub-function is the sum of the integral of the error modulus $e\left(t\right)$ multiplied by the time from the start of the transient process $t$ during the process $T$, and penalty additives $F_i\left(T\right)$. One of the most effective penalty additives is the positive part of the product of the error and its derivative [11,12,17,18].

The first member of the cost function has the following form:

$$F_1\left(T\right)={\int }_0^T\left|e\left(t\right)\right|tdt.$$

(6)

The second term is determined by the following relationship:

$$F_2\left(T\right)=k_w{\int }_0^T\left[\mathrm{m}\mathrm{a}\mathrm{x}\right\{0,e\left(t\right)\dot{e}\left(t\right)\left\}\right]dt.$$

(7)

Here

$$\max \left\{0,M\right\}=\left\{\begin{array}{c}0,ifM<0\\ M,ifM>0\end{array}\right..$$

(8)

$$\dot{e}\left(t\right)={\displaystyle \frac{e\left(t\right)}{dt}}.$$

(9)

Thus, the cost function has the following form:

$$F_{\mathsf{\Sigma }}\left(T\right)={\sum }_{m=1}^n{F}_{1m}\left(T\right)+F_{2m}\left(T\right).$$

(10)

The following modifications of the cost function components can also be used:

$$F_{1\mathrm{*}}\left(T\right)={\int }_0^T\left|e\left(t\right)\right|t^{2}dt.$$

(11)

$$F_{2\mathrm{*}}\left(T\right)=k_w{\int }_0^T\left[\sqrt{\mathrm{m}\mathrm{a}\mathrm{x}\{0,e\left(t\right)\dot{e}\left(t\right)\}}\right]dt.$$

(12)

$$F_{2\mathrm{*}\mathrm{*}}\left(T\right)=k_w{\int }_0^T\left[\sqrt[4]{\mathrm{m}\mathrm{a}\mathrm{x}\{0,e\left(t\right)\dot{e}\left(t\right)\}}\right]dt.$$

(13)

We will show in the examples below for what purpose it is possible to use a fractional power of the values included in the objective function, in particular, the square root of the signal or the fourth root. The rationale for using the square root is that the product of the error and its derivative is a quadratic form, which is proportional to the square of the error amplitude. This is not desirable, because we optimize linear system. The result should be invariant to the value of the test signal. So using the square root is a natural solution that has not been previously published yet. Modeling has shown that this improves the performance of the cost function.

The rationale for using the fourth root is based on two arguments. The first argument is that if introducing a square root improves the result, it is useful to try to apply a higher order root by induction as well. The second argument for this approach is that such a nonlinear function increases the contribution of small signals. This allows for a reduction in small signal deviations from the desired value when large deviations are already sufficiently suppressed.

During modeling, it is necessary to feed a test signal to the system model, which represents a task for the system, that is, a vector of input signals describing the desired values of the vector of output variables of the control object.

Since independent control of each output is required, the individual components of this vector must be substantially different functions and there must be no linear relationship between them. Traditionally, when optimizing a scalar linear system, a signal in the form of a single step jump is used [11,12,17,18]:

$$v\left(t\right)=\sigma \left(2∆t\right)=\left\{\begin{array}{c}0,t<2∆t,\\ 1,t\ge 2∆t\end{array}\right..$$

(14)

Here $∆t$ is the time discretization step. For this reason, the task vector is proposed $v\left(t\right)$ form in the form of stepwise jumps, shifted relative to each other by half the time of the transient process simulation [58]

$$V\left(t\right)=\left[\begin{array}{c}\sigma \left(2∆t\right)\\ \sigma \left(T/2\right)\end{array}\right].$$

(15)

In this case, instead of a multiplier proportional to the time from the start of the transient process, a sawtooth signal should be used, linearly increasing from the moment $t=0$ to the moment $t=T/2$, after which it resets to zero and again linearly increasing from zero at the same rate.

A delay in the occurrence of a step jump by at least the value equal to the doubled time sampling step $∆t$ is a mandatory condition that we have not encountered anywhere in the literature, but on which we have grounds to insist. If a jump at the input occurs without such a delay, then the simulation results may not coincide with the practical results, since in this case the system processes not the jump, but the initial state of the object, the results of signal differentiation in the system in this case will differ. A reliable check of the quality of the system is based on its operation when step jumps are supplied, and not on the type of processing of the error of the initial state. In the practice of modeling and optimization, we have had cases when the results of modeling in these two non-coinciding modes were noticeably different.

VisSim software (in any version) contain a mixture of functions in the time domain and functions presented as transfer functions, i.e. in the Laplace transform domain. There is no contradiction in this, since no error will arise. The mathematical notation of equations in such a mixed form would be erroneous, but the mnemonic designation of the blocks used in programming the project is quite understandable and for this reason acceptable. The VisSim software, designed for modeling and optimizing automatic control systems with feedback, uses various mathematical models. It contains models of linear elements that are described by transfer functions in the Laplace transform domain, as well as models of nonlinear connections that are described in the time function domain, and other elements. The software does not require uniformity of terminology, since it is not necessary to create full mathematical model of the system in the form of the one-type equations in the standard form for all of then form as correct mathematical description. The model is created from individual elements, between which, when this model is graphically specified, the appropriate connections are established from the inputs to the outputs; different ways of describing individual elements are not a problem.

5. Statement of the Problem

Example 1. Let us consider an object whose transfer function has the following form:

$$W\left(s\right)=\left[\begin{array}{cc}{\displaystyle \frac{2}{1+2s+s^2}}·e^{-0.5s}& {\displaystyle \frac{2}{1+3s+6s^2}}·e^{-0.8s}\\ {\displaystyle \frac{0.5}{1+4s+2s^2}}·e^{-0.7s}& {\displaystyle \frac{4}{1+4s+2.5s^2}}·e^{-0.3s}\end{array}\right].$$

(16)

We consider such objects for the following reason. A set of simple hypothetical models most accurately reflects all aspects of the possible problems considered in this article. Since not the only favorable type of models was used, but the inverse variant was also investigated, then the view to the chosen part of the problem can be treated as the researched one. More complex objects in the class of ones considered in this article are characterized by positive poles of some polynomials in the denominators of most specific elements of the transfer function. We do not consider such problems in this article. It is in our plans for the future.

If the polynomial of the transfer function denominators has no positive roots, then any example of the same order (an object of dimension $2\times 2$, polynomials in the denominator of dimension 2, polynomials of order zero in the numerator) will not be more difficult than the examples considered in the article. If the polynomial of the transfer functions has positive roots, then slightly different methods are necessary, including, for example, non-linear controller or some others.

6. Methods of Problem Solving

When modeling in the VisSim software, the object is formed from four models, which are shown in Figure 2. It should not be surprising that the coefficients of the free terms in the denominators of all transfer functions are equal to unit. This is the canonical form of these functions, since if this were not so, then the numerator and denominator of these transfer functions would be sufficient to divide by these coefficients to bring them to the canonical form. In this form, the transfer function at $s=0$ becomes equal to the coefficient in the numerator.

VisSim has three built-in numerical optimization algorithms. The Powell method, the Polak-Ribiere method, and the Fletcher-Reeves method. These methods are well known enough to be described in this paper, and, moreover, the possibility of using these methods is entirely due to the developers of this software; a description of these methods is beyond the scope of the authors of this paper. We have repeatedly found empirically that the Powell method is the best for this class of problems, since in many cases when one of the other two methods fails, or even both, the Powell method still succeeds, resulting in a set of optimal parameters, i.e., parameters that give a minimum value to the cost function. If any of the other methods works, then it gives exactly the same result as the Powell method. For this reason, there is no point in using the other two numerical optimization methods. And there is no point in describing the algorithm of the Powell method.

Powell’s method is based on the calculation of tangents and the assumption that any line passing through an extremum point intersects tangents to points of equal value of the objective function at the same angle. The other two methods are based on the calculation of the gradient of the objective function, which has a chance of unstable operation near the extremum point. At these points, the gradient is small and the influence of noise or random factors can lead to wandering around the extremum point. However, if the problems discussed in Appendix A do not arise, then all three built-in optimization methods can give the same optimization result.

7. Results

7.1. Two-Channel Object with a Majorizing Main Diagonal

Figure 3 shows the design for optimizing the controller for this case. In this figure, the individual PID controllers and the transfer function elements of the plant are combined into composite blocks called “PID” and “Compound”, respectively.

Figure 3, like the following figures, requires some preliminary comments. All these figures were obtained by screenshot from the VisSim software, so they have not been edited. Rectangles with an asterisk symbol inside denote signal multiplication devices. Rectangles with the inscription «pow» denote power raising units, unless otherwise specified in the text, this is the second power. Rectangles with the inscription «derivative» are input signal derivative calculators, rectangles with the symbols «1/ S» are integrators, a rectangle with a step or sawtooth signal drawn inside is a signal generator of this type, a rectangle with the inscription «abs» is an absolute value calculator of the input signal, a rectangle with the inscription «cost» is the input of the optimization device, where the calculated cost function is fed, which makes it possible to program the formula for its calculation.

The remaining rectangles are bus labels; they give the name to the buses or the quantities used. For each identical mark, it is considered that they are all connected to each other, and only one such block can be supplied with a signal; the program does not allow supplying another signal to another such mark, which is a natural limitation. In addition, the signal value of each mark can be used as a gain factor in another place. The blocks that set the gain factor have the form of a pentagon. In order to set the value of the variable gain factor in this way, it is enough to enter the name of the corresponding bus in the field for recording its value. Also, instead of the bus name, an ordinary number can be written in this field. Decimal numbers in this program are shown in a special format: if the number is less than zero, then the zero before the decimal point is not written, in negative numbers in this case there is a hyphen before the point, denoting a minus. The exponent is written after the number and separated by the letter “e”. Different colors of the lines correspond to the coloring of the input arrows. In particular, the output signal of the first channel is fed to the oscilloscope through the red arrow, so it is displayed as a red line. The output signal of the second channel is fed through the blue arrow and is displayed by a blue line. In all graphs here and below, the abscissa axis shows time in seconds, and the ordinate axis shows the output signal in conventional units.

The output values of any controlled object are physical values. These can be voltage, displacement, heating, rotation angle, etc. Each such value is measured by the corresponding sensor and then compared with the task. When modeling, it is not customary to write the physical units of these values, because control theory does not consider the absolute values of the output values, but only their increments in relation to the initial equilibrium state, and they are most often normalized, that is, divided by some standard value, so that only increments in relative units are displayed on the graphs. For linear systems, this is absolutely normal, since scaling the signal does not change the nature of the graphs. We do not consider nonlinear systems in this article.

The monitors on the center of Figure 3 show the obtained values of the corresponding coefficients of the matrix PID controller.

All the above mentioned approaches to optimization were applied here.

The software for modeling and optimization is VisSim 5.0e. This program is optimal for such tasks. Unlike MATLAB Simulink this software simply does not allow modeling of objects that cannot be implemented in practice. This software takes up little space on the disk, and the files it creates also take up little space. The main advantage is the precise reproduction of the algorithms by which real digital regulators operate.

The test signal, as indicated above, is single step jumps, shifted relative to each other by half the simulation time (15).

The simulation conditions are as follows: time discretization step $∆t=0.01sec$ and simulation time $T=80sec$. The integration method is the simple Euler method.

The objective function (10) with its constituent functions (6), (7), (8), (9).

The optimization method is Powell’s method.

The method of analyzing the result is a visual assessment of transient processes, an assessment of the roughness of the solution by rounding the obtained coefficients to 1%, and, if necessary, also changing the parameters of the object model by values $\pm 3\%$ or more.

The method of improving the result is to change the weighting coefficients in the cost function. To eliminate overshoot, non-monotonicity, oscillations, it is necessary to increase the weighting function for the term that depends on the positive part of the product of the error and its derivative; to reduce the static error, it is necessary to decrease this coefficient or increase the modeling time. If this does not help to obtain the desired result, then also increase the modeling time, modify the cost function.

In this result we do not see overshoot in the transient process of each channel, but we do see the response of each channel to a task step in the other channel, and this response is about 10% of the corresponding step.

Modeling and optimization showed that further changes in the cost function, including changes in the weight coefficients of the terms dependent on the product of the error and its derivative, did not significantly improve the transient process. Any specific improvement in one of the parameters of the transient process is achieved only at the expense of worsening the other parameter. Namely: it is possible to completely suppress this reaction in the second channel by increasing such a reaction in the first channel and decreasing the response speed of both channels. Or it is possible to increase the response speed of the second channel and almost completely suppress this reaction, but at the expense of increasing such a reaction in the first channel to 18%. It is also possible to significantly improve transient processes in the first channel, but at the expense of significantly worsening the response of which channel: in it, the advising response to a jump in the task of the first channel increases to 18%.

It should be noted that a relatively successful solution to the problem is achieved due to the fact that the transfer function has advantageous combinations of elements. This property is hidden in the fact that, firstly, the diagonal elements contain transfer functions in which the transfer coefficients are significantly greater than similar coefficients in other elements of the matrix transfer function of the object. Secondly, these same elements have a smaller delay value. Thirdly, the denominators of these same elements has smaller coefficients at the terms of the second and third order.

For educational purposes it is advisable to use objects like Example 1. Also, if in a practical task this is what a designer should strive for and sometimes this can be achieved by correctly choosing the numbering of the inputs in accordance with the numbering of the outputs.

Let us show that the indicated relations simplify the solution of the problem. To do this, we will zero out the elements in the regulator matrix that are not on the main diagonal. To do this, it is sufficient to zero out the corresponding coefficients, as shown in Figure 3, and remove from the project the extra blocks «Parameter Unknown». As can be seen from Figure 3, which shows the resulting project, control remains successful and even the deterioration in control quality was not too great.

However, this is the case, since in this case the simulation time had to be doubled, and, accordingly, the transient processes doubled. If in Figure 3 the processes when applying a jump to the first input last 10 seconds, and when applying a jump to the second input last 20 seconds, then in Figure 4 the processes when applying a jump to the first input last 20 seconds, and when applying a jump to the second channel last 30 seconds.

7.2. Two-channel object with a majorizing small diagonal

Example 2. Let us consider an object whose transfer function has the following form:

$$W\left(s\right)=\left[\begin{array}{cc}{\displaystyle \frac{0.5}{1+4s+2s^2}}·e^{-0.7s}& {\displaystyle \frac{4}{1+4s+2.5s^2}}·e^{-0.3s}\\ {\displaystyle \frac{2}{1+2s+s^2}}·e^{-0.5s}& {\displaystyle \frac{2}{1+3s+6s^2}}·e^{-0.8s}\end{array}\right].$$

(17)

This example is derived from the object of Example 1, in which the first column is replaced by the second, and the second column is replaced by the first. As a result, the transfer function of the object has become maximally unfavorable according to the criterion considered above. The elements on the main diagonal are characterized by a large delay, smaller transfer function coefficients, and large values of the coefficients in the denominator polynomials.

This example is one of the most unfavorable options for objects from this class.

Figure 5 shows the result of an attempt to optimize the diagonal regulator controller for the object of example 2. The attempt was unsuccessful, the system turned out to be unstable.

Figure 6 shows the project and the result of an attempt to calculate a full matrix PID controller for the same object from Example 2. A stable system is obtained. However, the quality of the obtained system is low. The overshoot in the first control channel reaches 98%. The overshoot in the second channel is 12%. The overshoot in the first channel with a jump in the second channel is 50%, the overshoot in the second channel with a jump in the first channel is 98%, as well as the reverse overshoot, which is about 50%.

In this case, the number of oscillations of all overshoots and with overshoot is at least 4 distinct oscillations. In addition, the duration of the process is very long, it is more than 60 seconds to the level of 5% and long residual relaxation processes of at least 40 seconds are also observed in the area of a small error. In addition, in the second channel, a second wave of response to a jump in the task of the first channel is observed, the peak of which is long and is between intervals of 20 seconds and 40 seconds.

From this example it should be concluded that with an undesirable combination of all parameters of the object model, multichannel control can be a rather complex task. However, this case is obviously bad; in practice, such cases should not arise, since the designer should pay attention to the fact that in the secondary diagonal all elements of the matrix transfer function have a significantly greater value in the entire frequency range for the reasons stated above. In this case, the designer should change the numbering of the outputs or the numbering of the inputs, which are the same in this case. To do this the output the first channel should feed the second input, and the output of the second channel output should feed the first input. This will result in the matrix transfer function of such an object taking the form given in Example 1, and this problem, as we have seen, is solved quite simply and successfully. This example shows that along with formal methods of designing controllers, a preliminary analysis of the properties of the matrix transfer function of the object should be used in order to identify the strongest channels of influence. In this case, the first input has the greatest influence on the second output, and the second input has the greatest influence on the first output, which dictates the choice of the rule for numbering outputs and inputs.

This approach is intuitive and should be used at an early stage of design, i.e., preliminary. For example, we want to move a tower crane load toward two directions and we have two adjustments, one for each arm. Let the left adjustment move the load strongly vertically, but weakly horizontally, and the right adjustment move strongly horizontally, but weakly affect vertical movement, then, when creating an automatic system, we should select the left adjustment to control vertical movement, and the right adjustment to control horizontal movement. Full autonomy will ensure that the left adjustment will affect only vertical movement, but will not affect horizontal movement, and the right adjustment, accordingly, will control only horizontal movement, but will not affect vertical movement. If the designer takes the opposite approach and tries to make the left adjustment control the horizontal movement and the right adjustment control the vertical movement, then the problem will still be solvable, but it will be similar to the one considered in Example 2, i.e. the least influential inputs are chosen to control these quantities, which complicates the controller design and can negatively affect the quality of control.

7.3. Two-channel object without an obvious majorizing diagonal

Example 3. Let us consider an object whose transfer function has the following form:

$$W\left(s\right)=\left[\begin{array}{cc}{\displaystyle \frac{2}{1+3s+6s^2}}·e^{-0.8s}& {\displaystyle \frac{4}{1+4s+2.5s^2}}·e^{-0.3s}\\ {\displaystyle \frac{0.5}{1+4s+2s^2}}·e^{-0.7s}& {\displaystyle \frac{2}{1+2s+s^2}}·e^{-0.5s}\end{array}\right].$$

(18)

This transfer function is also obtained from the transfer function of Example 1 with appropriate permutations. In this case, in the first row, the most significant element is the element that is not on the main diagonal, and in the second row, the most significant element is on the main diagonal.

Figure 7 shows the structure project and the result of optimization of the full PID controller for such an object from example 3. This system is stable, but the response of the second channel to the step jump of the first channel is too large. It reaches 75%, and its duration is 30 seconds. The response of the first channel to the jump in the second channel’s task is 14%, its duration to the level of 50% of the maximum value is less than 5 seconds, but the duration of its noticeable part is about 40 seconds. In this case, control for each channel is carried out without overshoot; the duration of the processes is from 40 seconds (for the second channel) to 60 seconds (for the first channel).

It is also advisable to try to control this object using a simplified diagonal controller. Figure 8 shows the design and result of controlling the object of Example 3 using a diagonal controller.

In this case, the result is noticeably worse, but the resulting system is stable. The response of the second channel to the step jump of the first reaches 98%, and its duration is 50 seconds. The response of the first channel to the jump in the second channel task is 13%, its duration is about 40 seconds. In this case, control for each channel is carried out without overshoot; the duration of the processes is 50 seconds for each channel.

Example 4. Let us consider an object whose transfer function has the following form:

$$W\left(s\right)=\left[\begin{array}{cc}{\displaystyle \frac{2}{1+3s+6s^2}}·e^{-0.8s}& {\displaystyle \frac{-0.5}{1+4s+2s^2}}·e^{-0.7s}\\ {\displaystyle \frac{-4}{1+4s+2.5s^2}}·e^{-0.3s}& {\displaystyle \frac{2}{1+2s+s^2}}·e^{-0.5s}\end{array}\right].$$

(18)

This transfer function is derived from the transfer function of Example 3, with the difference that the elements in the minor diagonal are negative. Figure 9 shows the design and optimization result of the diagonal controller for this case.

In terms of process quality indicators, there is no significant difference with the result presented in Figure 8: cross-responses and emissions are comparable in size and duration.

Figure 10 shows the design and optimization result of the full controller for the object from Example 4. When compared with the result shown in Figure 7, it has a slight advantage, namely: in Figure 6, the response at the output of the second channel with a jump in the reference at the input of the first channel is 75%, and in Figure 10, the corresponding response is 61%. In terms of duration, these responses are comparable, both of them are about 50 seconds. The other indicators are also comparable.

7.4. Three-channel object with slightly dominant main diagonal

Note that in the literature we have not found reliable reports on solving problems of three-channel objects control, which would provide mathematical models of the object with numerical values of the parameters, and the control result. In this case, we could take these initial models and calculate regulators for them, and compare the result with the result that would be presented in such an article. In fact, such a problem is solved in this article for the first time.

Example 5. Let us consider an object whose transfer function has the following form:

$$W\left(s\right)=\left[\begin{array}{ccc}{\displaystyle \frac{2}{1+2s+s^2}}·e^{-0.5s}& {\displaystyle \frac{1.5}{1+4s+2s^2}}·e^{-0.6s}& {\displaystyle \frac{1.2}{1+5s+2s^2}}·e^{-0.7s}\\ {\displaystyle \frac{1.5}{1+3s+6s^2}}·e^{-0.6s}& {\displaystyle \frac{1.8}{1+4s+2s^2}}·e^{-0.4s}& {\displaystyle \frac{1.4}{1+2s+3s^2}}·e^{-0.6s}\\ {\displaystyle \frac{1.2}{1+5s+2s^2}}·e^{-0.7s}& {\displaystyle \frac{1.4}{1+3s+2s^2}}·e^{-0.6s}& {\displaystyle \frac{1.7}{1+3s+s^2}}·e^{-0.5s}\end{array}\right].$$

(19)

As a test task, we will also use a vector consisting of time-shifted single jumps, but for a visual result, we will feed a jump with a negative value to the second channel:

$$V\left(t\right)=\left[\begin{array}{c}\sigma \left(2∆t\right)\\ -\sigma \left(T/3\right)\\ \sigma \left(2T/3\right)\end{array}\right].$$

(20)

Figure 11 shows the design and numerical optimization result of the controller for this example. The cost function calculator for this case is shown in Figure 12. The parasitic response of the second channel reaches 45% with jumps in the first channel. The parasitic response in the first channel is 20% with jumps in the second channel. The parasitic response in the third channel with a jump in the second channel is about 36%, with jumps in the first channel –40%. The overshoot in each channel is negligible.

For comparison, Figure 13 shows the design and optimization result of the diagonal regulator. It can be noted that the response in the first channel with a jump in the second channel increased by one and a half times, but the response in the same channel with a jump in the third channel decreased significantly. Also, the parasitic response in the third channel with a jump in the first channel decreased by two times. Thus, there is no reason to claim that in this case the diagonal regulator demonstrates a worse result than the full regulator.

Figure 14 shows the project and the optimization result for the case where the weight coefficient of the second term in the cost function is reduced to an insignificant value $k_w=0.1$. In this case, the performance of each channel has increased more than twofold, parasitic responses have also increased by 1.5 – 2 times, and emissions of up to 10% have appeared in the first and third channels.

It is also possible to use the square root in the second term and the weighting factor in this case $k_w=0.1$. The result for this case with a full controller is shown in Figure 15, and Figure 16 shows the result for this case with a diagonal controller. In terms of transient response quality, these two results differ only slightly from each other, and at the same time, they both differ quite significantly from the results obtained with the other cost function parameters shown in the previous figures.

Thus, this example shows that in the case where the matrix transfer function of the control object is characterized by a slight predominance of diagonal elements in absolute value over the remaining elements, the result expressed as transient processes for each channel separately depends to a much greater extent on the parameters of the objective function during optimization than on whether a diagonal controller or a full controller is used. The use of non-zero elements in the matrix PID controller does not allow for a radical improvement in the quality of the transient process, which is determined mainly by a compromise between speed and overshoot for all channels as a whole.

To obtain more complete information and substantiate the conclusions, we will calculate the regulator for a slightly modified object.

7.5. Three-channel object with a clearly dominant main diagonal

Example 6. We simply increase the coefficients at the main diagonal elements of the matrix transfer function of the object by approximately one and a half times. We obtain an object whose transfer function has the following form:

$$W\left(s\right)=\left[\begin{array}{ccc}{\displaystyle \frac{3}{1+2s+s^2}}·e^{-0.5s}& {\displaystyle \frac{1.5}{1+4s+2s^2}}·e^{-0.6s}& {\displaystyle \frac{1.2}{1+5s+2s^2}}·e^{-0.7s}\\ {\displaystyle \frac{1.5}{1+3s+6s^2}}·e^{-0.6s}& {\displaystyle \frac{2.7}{1+4s+2s^2}}·e^{-0.4s}& {\displaystyle \frac{1.4}{1+2s+3s^2}}·e^{-0.6s}\\ {\displaystyle \frac{1.2}{1+5s+2s^2}}·e^{-0.7s}& {\displaystyle \frac{1.4}{1+3s+2s^2}}·e^{-0.6s}& {\displaystyle \frac{2.55}{1+3s+s^2}}·e^{-0.5s}\end{array}\right].$$

(21)

The cost function and test signal in this case are the same as in example 5.

Figure 17 shows the design and optimization result of the diagonal controller for example 6. The quality of the resulting transient process in all three channels has improved significantly in this case. The parasitic response in the first channel varies from 10% to 13%, in the second channel it remains within 20%, in the third channel it varies from 18% to 30%. There is no overshoot in any channel.

Figure 18 shows the design and optimization result of the full controller for example 6. The result is slightly better than in Figure 17. All parasitic responses have become less than 10%, and the responses of the first and second channels to the task jump in the third channel are even less than 5%, and there is no overshoot in each channel either.

Based on this example, it can be argued that if in the matrix transfer function of the control object there is a significant (2 times or more) predominance of the values of the diagonal elements over other elements, then with a diagonal controller, control of acceptable quality is achieved, and with a full controller this quality can be noticeably (twice) improved.

The indicated predominance can be called a majorizing property, which consists in the fact that each diagonal element is significantly larger in absolute value than the other elements in the same row of a given matrix.

7.6. Three-channel object without majorizing main diagonal

To make sure that majorization is only important within each row, consider another example based on the transfer function from Example 6, namely, we decrease the value of all coefficients in the first row by 2 times and increase the value of all coefficients in the second row by 2 times.

Example 7. We obtain an object whose transfer function has the following form:

$$W\left(s\right)=\left[\begin{array}{ccc}{\displaystyle \frac{1.5}{1+2s+s^2}}·e^{-0.5s}& {\displaystyle \frac{0.75}{1+4s+2s^2}}·e^{-0.6s}& {\displaystyle \frac{0.6}{1+5s+2s^2}}·e^{-0.7s}\\ {\displaystyle \frac{3}{1+3s+6s^2}}·e^{-0.6s}& {\displaystyle \frac{5.4}{1+4s+2s^2}}·e^{-0.4s}& {\displaystyle \frac{2.8}{1+2s+3s^2}}·e^{-0.6s}\\ {\displaystyle \frac{1.2}{1+5s+2s^2}}·e^{-0.7s}& {\displaystyle \frac{1.4}{1+3s+2s^2}}·e^{-0.6s}& {\displaystyle \frac{2.55}{1+3s+s^2}}·e^{-0.5s}\end{array}\right].$$

(22)

Figure 19 shows the design for optimization and the result for this case. The result is slightly worse in quality than the result shown in Figure 17. This is probably due to the change in the relative influence of each channel on the other channels, as well as the effect of the result on the cost function.

If designer is not satisfied with the quality of the transient process in a specific channel, then an additional weighting factor can be introduced into the corresponding term in the cost function depending on the error in this channel. For example, in this case, the response quality of the second channel is noticeably worse than the response quality of the other channels. The best response is in the first channel. Designer can introduce a factor of 3 into the second channel, a factor of 2 into the third channel, and leave a factor of 1 in the first channel.

Figure 20 shows the result. As we can see, we managed to reduce the parasitic responses in the third channel by increasing the parasitic responses in the first channel. This further demonstrates the dominant influence of the cost function during optimization and its numerical coefficients on the result.

8. Ensuring the requirements of completeness of the influence of test tasks

If we analyze the set of test signals used in the optimization, we can see that they are somewhat incomplete. Note that with a step jump at one of the inputs, there is no jump at the other input. The optimization procedure provides better transient processes in these situations. However, the possibility of a simultaneous jump of the reference at both inputs, which can act together in such a way that an overshoot of an inadmissibly large value will occur, is not considered. Such a situation can occur both when the reference jump signs coincide and in the case of opposite signs. To ensure complete completeness of all possible signal combinations, it is proposed to include in the set of test signals also the simultaneous supply of in-phase jumps to all inputs and the simultaneous supply of antiphase jumps to all inputs. For an object with a dimension, $2\times 2$this gives two more jumps, and for an object with a dimension, $3\times 3$this implies many more jumps in various combinations: in-phase and out-of-phase jumps on each two inputs in pairs, for a total of 6 jumps, also in-phase jumps on all three inputs, and in-phase on two inputs with out-of-phase on the third input, for a total of four options. Thus, this gives a total of 13 different options. We will demonstrate this method using the object with a dimension $2\times 2$from Example 1 and from Example 2.

Example 8. Let’s consider the object from example 1 and apply the following set of test tasks to it:

$$V\left(t\right)=\left[\begin{array}{c}\sigma \left(2∆t\right)-\sigma \left({\displaystyle \frac{T}{2}}\right)+\sigma \left(3T/4\right)\\ \sigma \left(T/4\right)-\sigma \left({\displaystyle \frac{T}{2}}\right)-\sigma \left(3T/4\right)\end{array}\right].$$

(23)

Figure 21 shows the design and optimization result for this case using the simplest cost function consisting of the integral of the sum of the squares of the product of the error modulus and the time since the start of the transient process.

At the same time, two large parasitic jumps in the second channel remained, namely, when the task jumps only on the first channel and when the task jumps unidirectionally on both channels at the same time. But the other components of the complex transient process improved. It is important that in this case there is no overshoot in the second channel, and the parasitic jump during task jumps in the first channel is also not large, only 10%. For comparison, we will show the result when only one jump is used on each channel. This result is shown in Figure 22. In this case, a similar jump in the second channel is 25%. From this we can conclude that considering the set of all four types of jumps is more representative and in this case, formally, the same cost function works differently.

It is noteworthy that at the moment in time $t=80sec$ in Figure 21 the overshoot in the first channel is 60%, from which it follows that the simultaneous synchronous jump in the task on both channels is a strong exciting factor for the first channel. At the moment of time $t=120sec$ the responce at the outputs of both channels are insignificant, which indicates that the antiphase jumps of the steps are not an exciting factor for both channels.

Finally, it is important to use both methods of improving the transient process at once, namely, the test signal is as in Example 8, and the cost function is full, as in Figure 9.

Example 9. Figure 23 shows the design and the result of optimization using four jumps in each channel and a complex cost function that helps reduce overshoot. The result can be called ideal: there is no overshoot in all cases for all channels, the parasitic response in both channels is negligible. Thus, complete consideration of all variants of the system input signals significantly increases the efficiency of the controller optimization design. It is enough to compare this result with the result shown for the same case in Figure 3 and Figure 4 - there the parasitic jumps in each channel were 10% with a duration of about 10 seconds.

For this reason, it is worth considering how the controller optimization would be performed in this case for the object obtained from the object of Example 1 by replacing the first output with the second, and the second output with the first. In this case, considered in Example 2, the result shown in Figure 5 is extremely poor, with a large overshoot and many oscillations in each channel.

Example 10. Consider the object of example 2 with test signals as in example 8.

Figure 24 shows the design for this case and the result of numerical optimization. The result exceeded all expectations. Despite the fact that in this case the elements of the main diagonal in the entire frequency range are smaller in absolute value than the absolute value of the elements of the non-main diagonal, the optimization result shows excellent quality of transient processes: there are no spikes, and the parasitic responses in each channel with jumps in the other channel are also negligibly small in value, no more than 1–2%.

This example particularly clearly shows that when completeness requirements are presented for the quality of transient processes, including the requirement that not only transient processes for each channel separately in the absence of jumps in the task for other channels have a sufficiently high quality, but also that processes with simultaneous in-phase and simultaneous antiphase jumps at all inputs of the system, transient processes would also be characterized by high quality: little overshoot or its complete absence.

9. Discussion

As a result of modeling and numerical optimization of two-channel full and incomplete PID controllers, the following conclusions were made:

Particular attention should be paid to the ratio of the values of the main diagonal elements of the object’s transfer function matrix to the remaining elements of this transfer function. Signs of a larger value are: a) a larger gain value; b) a smaller delay value; c) a smaller value time constants in the polynomial in the denominator. If these requirements are not met even partially, it is advisable to see whether they will be met to a greater extent if the numbering of the outputs is replaced by a different correspondence between the outputs and inputs.

When optimizing, it is recommended to use step jumps with a time shift by the corresponding fraction of the simulation time as a test task. In particular, for a dimension object, $2\times 2,$it is recommended to shift the jump on one of the channels by half the simulation time.

The very first jump must be shifted in time by an amount no less than twice the duration of the time quantization step.

The approach to solving the optimization problem proposed in this article recommends use the complex cost function (10) as the objective function.

If optimization does not ensure that transients are completed, it is recommended to try increasing the simulation time.

The use of a more complex test signal in the form of a sequence of single step jumps of different polarity, which contains separate jumps to each input, an in-phase jump to both inputs, and antiphase jumps to both inputs, is proposed and substantiated. It is shown that the use of such a test signal, while observing all other optimization conditions, can significantly improve the quality of transient processes, which in our two examples can be called ideal. Such a test signal with its experimental substantiation is proposed for the first time.

The simulation confirmed that if the matrix transfer function of the object is characterized by the predominant absolute value of the elements on the main diagonal, then even a diagonal controller allows for quite acceptable autonomy of the channels, i.e. separate control of the output values of each channel of the controlled object. In this case, a complete controller can give a better result. However, this improvement is qualitative, not fundamental. The properties of the matrix transfer function may be unfavorable, which is expressed in the fact that the elements of its main diagonal are commensurate in value with the other elements in the entire frequency range. Then the quality of the result of solving the optimization problem will be significantly worse, since control using a complete controller can be made stable, but this can lead to large overshoot values, large channel output responses with a step change in the task in other channels, as well as a large number of oscillations in both types of these transient processes. In some cases, an incomplete controller can provide stable control.

Apparently, with further growth of the number of channels, the design of a complete regulator becomes not only more complex, but also less and less expedient, whereas in the case of ensuring a given property of the matrix transfer function of the object, a diagonal regulator will be sufficient. In this case, the complexity of the design problem diagonal controller will not increase so sharply with the growth of the number of channels, since the problem actually degenerates into the problem of designing scalar controllers, the number of which corresponds to the number of channels in the object. It is possible that for this reason there are practically no publications describing the results of designing controllers for multi-channel objects with more than three channels, and the number of articles on the design of three-channel controllers is negligibly small.

10. Conclusion

The research results presented in this article offer designers of multichannel systems a complete methodology for designing sequential PID controllers, which, if necessary, can be applied to more complex controllers, including those with additional channels of double integration and (or) additional derivation. More complex controllers are hardly justified for solving such problems. The methodology consists of choosing a software optimization tool, an optimization method (it is recommended Powell’s method), the integration method (the simple Euler method is recommended), the choice of the objective function (recommended according to the relation (10) of two terms), a set of test signals of the problem in the form of single jumps shifted by equal time intervals and their superposition, providing: a) separate jumps at each input in the absence of changes at other inputs; b) simultaneous in-phase jumps; c) simultaneous antiphase jumps. The technique is also supplemented by an objective function based on the integral of the error modulus multiplied by a sawtooth signal with the addition of the integral of the square root of positive part of the product of the error by its derivative. The method of increasing the influence of one of the two specified terms in function (10) is the selection of the best weight function based on the analysis of the optimization result. Also, in case of failure, it is recommended to increase the simulation time.

Of course, any method has not only advantages, but also its limitations. The limitations of any method are revealed when solving specific problems. The numerical optimization method has fewer limitations than other methods for solving problems of this class. Unlike the tabular methods and the analytical methods, the numerical optimization method not only easily allows for delays and even nonlinear elements in the object model, but is also not much more complicated in these specific cases. The only difference in the presence of nonlinearity is that in this case, using as a task single step jumps is unjustified. In the presence of nonlinearity, the test signals should, if possible, coincide with the actually expected signals in the real system. The advantage of the numerical optimization method over engineering methods is its lower labor intensity due to the transfer of a large number of calculations to the computer, while the developer does not have to calculate practically anything. The method will not lead to the expected results if they are in principle unachievable for an object with a given mathematical model and given requirements for quality indicators. Unlike the Monte Carlo method used in the MATLAB Simulink application, the numerical optimization method does not find a solution that fits into the specified boundaries, but the best solution from the standpoint of the specified cost function. If the boundaries are set too broadly, the solution will be significantly worse than the best, if they are set too narrowly, the solution may be absent. Therefore, when using the Monte Carlo method, a method for determining the capabilities of the system for a given object is required; in the numerical optimization method for a given cost function, this is not required; the result is always optimal within the given modeling conditions and values of the weighting coefficients.

Financing: The work was carried out with the support of the Ministry of Education and Science of the Russian Federation (within the framework of state assignment No. 075-00682-24) and using data obtained at the unique scientific installation “Seismic and infrasound complex for monitoring the Arctic permafrost zone and complex for continuous seismic monitoring of the territory of the Russian Federation and adjacent territories of the world”.