Table of Contents Author Guidelines Submit a Manuscript
Journal of Control Science and Engineering
Volume 2008 (2008), Article ID 723292, 7 pages
http://dx.doi.org/10.1155/2008/723292
Research Article

Actuator Fault Diagnosis with Robustness to Sensor Distortion

INRIA-IRISA, Campus de Beaulieu, Rennes Cedex 35042, France

Received 1 April 2007; Revised 15 August 2007; Accepted 11 November 2007

Academic Editor: Jakob Stoustrup

Copyright © 2008 Qinghua Zhang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Actuator fault diagnosis is often studied under strong assumptions on available sensors. Typically, it is assumed that the sensors are either fault free or sufficiently redundant. The purpose of this paper is to present a new method for actuator fault diagnosis which is robust to sensor distortion. It does not require sensor redundancy to compensate sensor distortion. The essential assumption is that sensor distortions are strictly monotonous. Despite the nonlinear and unknown nature of distortions, such sensors still provide useful information for fault diagnosis. The robustness of the presented diagnosis method is analyzed, as well as its ability to detect actuator faults. A numerical example is provided to illustrate its efficiency.

1. Introduction

In the design of modern control systems, fault diagnosis is often considered for component faults, sensor faults, and actuator faults. Various methods for fault diagnosis are generally based on the processing of sensor signals [13]. Many methods for actuator fault diagnosis assume reliable fault-free sensors, see, for instance, [4, 5]. For methods simultaneously dealing with actuator and sensor faults, it is typically assumed that there is sufficient redundancy among the sensors such that at any moment, the valid sensors can provide the information required for faults diagnosis. For mass production, it is important to use as less sensors as possible to reduce production cost. In such situations, fault diagnosis cannot uniquely rely on redundant sensors. Moreover, still for the purpose of cost reduction, sensors used in mass production may be highly nonlinear. If the sensor nonlinearity is well known, it can often be electronically compensated. Unfortunately, sometimes sensor nonlinearity varies within the production, and for each piece in use, the nonlinearity varies during its normal life duration. For example, most oxygen sensors used in cars equipped with catalytic converters are highly nonlinear. They roughly indicate if the oxygen concentration is over or below a reference value. Nevertheless, such sensors are sufficient for the purpose of engine control. It is then important to develop fault diagnosis methods relying on the same sensors. In this paper, unknown nonlinear behaviors of sensors are generally called sensor distortion.

The purpose of this paper is to present a method for actuator fault diagnosis which is robust to sensor distortion. It is assumed that each sensor can be affected by an unknown and arbitrary, but strictly monotonous, nonlinearity. The monotonousness is a weak assumption since any non monotonous distortion would make the sensor information useless. Remark that saturation is not a strictly monotonous nonlinearity, and thus is excluded in the proposed method. Nevertheless, a correctly working sensor should not be saturated when it is in its normal working range.

Robust methods for fault diagnosis has been studied in different contexts by many authors. In [6], a method based on adaptive wavelet analysis is proposed. The method studied in [7] uses generalized frequency response functions. Kullback discrimination information is used as a fault detection index in [8]. In [9], some faults of induction motors are detected by Fourier analysis of frequency signatures. (See the references cited in these publications for more information.) Compared to these existing results, the novelty of the method presented in this paper is its ability to deal with unknown nonlinear sensor distortions without requiring sensor redundancy.

The preliminary results presented in [10] are further developed in this paper, in particular, the analysis of the proposed robust detection method is completed with the characterization of the faults which can be detected.

The paper is organized as follows. The problem considered in this paper is formulated in Section 2. Fault detection and isolation are, respectively, studied in Sections 3 and 4. A numerical example is presented in Section 5. Some concluding remarks are given in Section 6.

2. Problem Statement

The considered linear state-space system with sensor distortion is formulated as ?????????????????????????????????????????????=????+??diag??+??,(1a)??(??=????+??,(1b)=h,(1c)??(??)?R????(??)?R?? where ??(??)?R?? is the state, ??(??)?R?? the input, ??(??)?R?? the output before sensor distortion, ??(??)?R?? the output after sensor distortion, and diag(??(??)) and ??(??) represent bounded modeling uncertainties. The notation ???R?? denotes the diagonal matrix formed by the components of the input vector ??, and the vector ??0 is a coefficient vector introduced to describe the efficiency loss of actuators (multiplicative actuator faults). For fault-free actuators, (??) takes the nominal value h.

For notation simplicity, the parenthesis ???? of the time-dependent variables will not be written unless necessary.

Sensor distortion is modeled by the component-wisely defined nonlinear function ????: let ?? and ?? be, respectively, the ??th component of ??=1,,?? and ????=h????????.(2), h??:R?R, then ??1

Assumption 1. Each sensor distortion ??2?R is an unknown, but a strictly monotonously increasing, function. In other words, for any ??1<??2?h?????1?<h?????2?.(3), ??, ??

Assumption 2. The system dynamics matrix ??????=h+??,(4) is asymptotically stable, that is, the eigenvalues of ?? have negative real parts.

Remark 1. No additive uncertainty is assumed in the sensor distortion equation (1c) since the arbitrary unknown monotonous nonlinear function can take into account some sensor distortion uncertainty. On the other hand, if it is assumed that ??=h-1?h???????+??+??-????,??=????+??,(5) where ????=h??(6) is some additive uncertainty, define??then (4) can be replaced by hwhich is in the form of (1c). Of course, in order to limit ??, some regularity of h?? should be assumed.

With the above formulation, the problem considered in this paper is the detection and isolation of multiplicative actuator faults, modeled as changes in the coefficient vector ????=h??(????), despite the unknown sensor distortions.

3. Fault Detection

The main difficulty of the problem formulated in the previous section is caused by the unknown sensor distortions. The key question is how to use the information provided by such sensors. Because of the unknown nature of ????, for each measured value of ????, the corresponding value of h?? is completely unknown, and even the sign of ?? is unknown. The strict monotonousness of ?? assumed in Assumption 1 is not helpful in this aspect. However, it is important to make the following observation. For any two time instants ???sign??????-????????????=sign??????-?????????.(7) and ????, Assumption 1 implies that ????In other words, the relative sign of ?? at different time instants is known from the sensor output ?? measured at these time instants. Remark that |????(??)-????(??)| and ????(??)-????(??) are two arbitrary and independent time instants, and either one can be earlier than the other one. Since the absolute value ??(??) is completely unknown, the relative sign is thus the only information about ??(??) provided by the sensor output. This information will be the basis for the design of fault detection and isolation algorithms in this paper. Such information can also be used for the identification of Wiener systems [11].

Let ??(??), ??(??), ??(??), and ??(??) be, respectively, the Laplace transforms of ??(??), ??(??), ????????=??????-??-1????????????diag??+??????-??-1??????????+??.(8), and F????=??-1??[??????-??-1????????????????diag],(9)=??-1??[??????-??-1??????????+??],(10). It is then derived from (1a), (1b), and (1c) (by assuming zero initial state) that ??-1 Define ????????????????????????=F??+??(??),(11a)=h.(11b)F(??)?R??×?? where ??(??) is the inverse Laplace transform operator. Then ??(??)?R????(??)

Notice that ??(??) depends on ?? and can be computed through (9), whereas F depends on modeling uncertainties F, ???? and is unknown.

Remark 2. Zero initial condition of the state ?? has been assumed when (8) was derived. The asymptotic stability condition (Assumption 2) is not explicitly used in the above reasoning, but it ensures well-behaved computation of F through (9). To some extent, this stability condition allows also to tolerate nonzero initial states which are asymptotically forgotten. If the system was not asymptotically stable, then, in principle, an observer should be used in the computation of ????. However, it is difficult to design observers with sensors distorted by unknown nonlinear functions.

Let ?? be the ??th row of ?? and ????????=??????????+????????.(12) be the ??th component of ??, then the ???sign??????-??????????????????-?????????=0.(13)th row of (11a) writes????(??) For any two time instants ????(??) and ??, it is obvious that ???sign??????-???????????????????-???????????+????????-?????????=0.(14) Substitute ???sign??????-???????????????????-???????????+????????-?????????=0,(15) and ???-sign??????-??????????????????-??????????????=sign??????-??????????????????-?????????.(16) with the last equality, then, for constant ??0, ?? Remind the relative sign equality (7), then ????(??,??) or equivalently???????????,??=-sign??????-??????????????????-???????????0.(17)

Let ???????????,??=sign??????-??????????????????-?????????.(18) be the nominal value of |????(??)|=?? (corresponding to fault-free actuators), a residual ?? for actuator fault detection can be generated as ?? It then follows from (16) that, in the fault-free case, ????(??,??) This result leads to the following proposition.

Proposition 1. If ??=??0 for some constant ????????,??=2??.(19) and any ????(??)?????(??), then the residual ??????=|||????,??????|||=|||????,????????-????????|||=|||????????|||+|||????????|||=2??.(20) defined in (17) satisfies, in the fault-free case, that is, ????(??)=????(??), the inequality ????????,??=0=2??.(21)

Proof. Let us first consider the case ????(??), then it follows from the inequality (18) that ??(??) Now for the case ??(??), it follows trivially from the inequality (18) that ??(??)

Remind that ??(??) is a component of the variable ??(??) depending on the bounded modeling uncertainties ??(??) and ??(??), as defined in (10). The boundedness of ??(??) and ????(??,??), together with Assumption 2, implies the boundedness of 2??. For practical convenience, the residual threshold should be directly derived from the assumed bounds of (??,??) and ???0. It would require nontrivial error bound propagation through (10). Techniques of interval analysis or set-membership computation [12] can be applied for this purpose. This topic is not further discussed in this paper.

Proposition 1 guarantees the absence of false detection if fault detection is made by comparing the residual ??0?0 with the threshold ?????0. However, it does not tell what are the faults which can be detected with such a decision rule. In general, robust fault detection methods are based on conservative decision rules, preventing the detection of some faults. In publications about robust fault detection, typically robustness results are provided, but the analysis about the set of faults which can be detected is usually absent, because it is often difficult to characterize the detectable faults in a robust detection framework. In contrast, for the method proposed in this paper, the faults which can be detected are clearly characterized as follows.

Proposition 2. Assume that the system matrix pair ?? is controllable, ??=????0, and ??(??). For the faulty actuator parameter vector ??, if there does not exist any positive real number ?? such that ????????,??>2??(22), then there exist an input signal ??=1,,?? and time instants ?? and |????(??)|=?? such that ????(??,??) for each ??=1,,??, where ????????=??(-1??????-0??-1??0),(23) is a positive constant such that v????=??????.

Proof. This proof applies to the residual ?? for each ??????????-??????=??(0??-1????0????????????)=??(1-0??-1??????-1????0??)=0,(24).
Define the vector |????0??|=???0????? where the norm |????0??????|=0????????(25) and ??=????0 is a positive number to be specified later. Then ??>0 where the last inequality follows from the fact that ??????????-??????=??(0??-1????0??)>0.(26).
Notice that the equality ??????0??????-????=??(0??-1????0??)<0.(27) would imply sign[??????]??????0<0.(28) for some ??, that is, the case excluded in Proposition 2. Therefore, ??0
Similarly, it is also shown that ?? It then follows that ???? For given values of ?? and F, the value of ??(??) in (23) is chosen large enough such that ????(??)-????(??)=????
Now, let us consider ??, the ??th row of the matrix ??sign????????-?????????????????????-???????????0<0.(30) defined in (9). If there is an input signal (??,??) such that |||(????(??)-????(??))??|=|??????|>2??=|????(??)-????(??)| for two time instants (????(??)-????(??))??+(????(??)-????(??)) and (????(??)-????(??))??, then the inequality (28) leads to ??sign????????-??????????????+??????-????????????=sign????????-????????????.(31)
Notice that the existence of such input signals is ensured by the controllability of the matrix pair ??sign????????-??????????????+??????-???????????????????-???????????0<0,(32).
Because |[????(??)-????(??)]??0|=|??????0|>2??, the sign of ???sign??????-??????????????????-???????????0<-2??.(34) is determined by the sign of ????????,??>2??.(35), thus ?? This last equality, together with the inequality (30), leads to ????(??,??) or equivalently ?? Remind that ????(??,??), then ?? Therefore, the residual as defined in (17) satisfies F(??)

If ??(??) is the current time instant, then 1,2,,?? can be computed for different past time instant 0<??<??. In order to reduce the effect of modeling uncertainties, ??=1,2,,?? can be averaged over different values of ??=??+1 in a sliding window. The computation of the residul from sampled signals is summarized as follows.

Residual Generation Algorithm Summary
Assume that ???????????,??=-sign??????-??????????????????-???????????0,????????=1????-1???=??-?????max?????.??,??,0(36) and ????(??,??) are sampled1 at discrete time instants ??. Choose the sliding window length ??0 for residual averaging. For each sensor number h, a residual is computed, for ??, through the formulas ???R Notice that the max function is used to exclude negative values of ??1=????0 in the computation of the average value.

4. Fault Isolation

After the detection of an actuator fault, the purpose of fault isolation is to figure out which actuators are faulty. In terms of (11a) and (11b), it amounts to deciding which components of ??0 have deviated from the nominal value F(??).

It should be first remarked that, because of the arbitrary unknown function ??(??), it is not possible to detect or isolate any proportional changes in all the components of ??. In other words, for any value ??, the parameter vector ??×?? cannot be distinguished from ???? based on the known signals ?? and ??. For the same reason, if all the components of ????, except one, have changed, it is not possible to determine which one has not changed.

After having clarified the limitation related to unknown sensor distortion, let us look for an algorithm for fault isolation. The basic idea is to design residuals similar to (17), but capable of rejecting some actuator faults. Each designed residual should be insensitive to some of the possible actuator faults, whereas sensitive to the others. The actually occurred fault can then be isolated by comparing such residuals.

Let ???? be a permutation matrix, that is, a matrix obtained by permuting the rows of the ???? identity matrix ????=??-???? (remind that ?? is the number of actuators). Divide ??????????+??????????=????.(37) into two sub-matrices ?????? and ??????, respectively, composed of ?? and ????????????=F????????????????+F??????????????+??.(38) rows of ????. Then it can be easily verified that ??

The notations ??????=??????0(39) and ???????????,(40) will be used to select, respectively, the assumed faulty components and sound (or fault-free) components of ????????????=F??????????????+F??????????0????+??.(41). It is derived from (11a) that ?? If (????,????) is assumed to select the components of ????=?????? corresponding to sound actuators, then, even after the occurrence of actuator faults, the equality ???? still holds. Define ???? then ??

For the purpose of fault isolation, different partitions of ?? into ???sign??????-???????????????????-???????????????????+?????????-?????????×????????????0+????????-?????????=0.(42) should be considered. For each particular partition, a residual will be designed to be insensitive to changes in the corresponding subvector F(??). Such a residual is said to be rejecting changes in ??(??). The rejection method used in the following is the estimation of 1,2,,?? from measure signals.

For two time instants ????, ????(??)-????(??), Assumption 1 implies (by rewriting inequality (15)) that?? For a given set of signals ??, 1,2,,?? sampled at discrete time instants ????(??)<????(??), the value of ????(??)<????(??) can be estimated by minimizing the error term ????(??)<????(??) in the inequality (42) where the time instants 1,2,,??, ???? are replaced by all pairs among the sampling instants ????(??).

Because ?? and (????,????) imply ????=??????, there would be too much redundancy if all the possible pairs among ??=1,2,,?? were considered in (42). In order to reduce redundancy, for each sensor ????(????1)=????(????2)=?=????(??????).(43), the data samples are sorted according to the values of min????max1=??=??,1=??=??-1????????,(44), and only the neighboring pairs are considered. The algorithm is summarized as follows.

Residual Generation Algorithm Summary
For each chosen partition of ??=1,2,,?? into ??=1,2,,??-1, the residuals rejecting changes in [????(??????)-????(??????+1)]??????????+[????(??????)-????(??????+1)]????????????0+????????=0.(45) are computed as follows. For each sensor number ?????1?,?????2???,,??????-1,??=1,2,,??,(46), sort the data such that ???? Solve the constrained optimization problem sign[????(??????)-????(??????+1)] subject to the constraints, for ?? and ??????=??????0, |????(??)|=?? The corresponding sequences ?? are the residuals rejecting changes in ??.

Remark that the inequality (45) corresponds to the inequality (42) where the omitted ???????? is always negative due to the sorted sequence, and accordingly, the inequality changes from “=” to “=.”

The constrained minimization (44)-(45) can be easily reformulated in the form of a standard linear programing problem. There are efficient numerical algorithms for its solution [13].

Proposition 3. If the true parameter vector ????????=2??.(47) governing (11a) satisfies [????(??????)-????(??????+1)]????????????+[????(??????)-????(??????+1)]????????????+????(??????)-????(??????+1)=????(??????)-????(??????+1)=0,(48), and if ?????? for some constant ??????+1 and any ????(??), then the residual ??????=??????0, solution of the constrained optimization problem (44)-(45), satisfies ????

Proof. The proof of this result is quite straightforward. Let us first derive from (11a), (11b), and (37) that ?????? where the indices ????????, ????(??????)-????(??????+1) come from the sequence sorting ???? as in (43).
It is assumed that ????????. Then the above inequality shows that there exist a value of ???????? (equal to ????(??????)-????(??????+1)) and a value of ????????=????(??????)-????(??????+1)=2??.(49) (equal to ??????=??????0) such that the inequality (45) is satisfied. These values of ??????=????????0 and ???R do not necessarily correspond to the solution of the constrained optimization problem (44)-(45), however, the optimal solution certainly has a value of ???? not larger than ????. Therefore, ????

Property (47) has been proved under the assumption ??????=????????0. It can be shown that (47) holds also if ?? for any ???? because of the unknown sensor distortion function.

For fault isolation, various matrices ?? are assumed and the corresponding residuals are computed. Then property (47) can be used to decide if each assumed ??????=??????0,(50) is correct or not. Alternatively, for different matrices ?????????? of the same size, the values of the residuals can also be compared for fault isolation.

Keep in mind that fault isolation cannot distinguish the cases such that ?????. For reliable fault isolation, the following assumption is required.

Assumption 3. There exists a permutation matrix ???? such that, for ??????????=??????0,(51) composed of some rows of ??, ??????=4.38??-20????33??+11????+4.42??-22????44??+12????+7.20??19??+13????.(52) and there does not exist any ° (except for h????=1001+??-0.2???-320?+10(53) composed of a subset of permutated rows of ??(??)) such that ° for some scalar real value ??=100.

5. Numerical Example

In this section, the presented fault diagnosis method will be illustrated with a simulated distillation column.

In [14, page 223], a distillation column is modeled with a transfer function matrix with three control inputs: top draw flow rate, side draw flow rate, and bottom temperature control, and three outputs: top draw composition, side draw composition, and bottom reflux temperature. In order to illustrate the possibility of fault diagnosis with only one output sensor, let us consider the model relating the three inputs and one of the outputs, the bottom reflux temperature. The transfer function model is Note that the known time delays at the inputs do not cause any serious difficulty for fault diagnosis, though the on-line computation has to be delayed accordingly.

Numerical simulation is first made in continuous time. The three input variables are randomly drawn with uniform distributions ranged within the intervals [20, 40]?(mol/min), [10, 30]?(mol/min) and [10, 20] (C). The simulated noise-free output is disturbed by an additive bandlimited noise with noise power 0.01 and unitary sample time. Then the sensor distortion function is applied. In this monotonously increasing function, the additive constant 10 has been added to illustrate the robustness of the proposed method to sensor bias (shifting error). Finally, the distorted sensor output sampled at the period of one minute is corrupted by a random noise uniformly distributed within the interval [-0.05, 0.05] (C). The simulated data is recorded during 2000 minutes after an initial simulation of 1000 minutes to avoid the initial transient period. At the beginning, all the actuators are fault free. At the 1000th minute (of the recorded duration), a factor of 0.8 is applied to the first actuator, simulating an actuator fault.

The residual for fault detection, computed with the averaging window length , is illustrated in Figure 1. In Figure 1(a) showing the residual in full scale, there is a strong transient behavior of the residual after occurrence of the actuator fault (at the 1000th minute). This transient behavior is in favor of a fast detection of the actuator fault. In order to view better the residual outside the transient period, it is plotted in Figure 1(b) in a finer scale. The detected fault is clearly confirmed by the residual after the transient period.

fig1
Figure 1: Fault detection residual. The same residual is plotted twice at different scales. The time unit is the minute.

For fault isolation, three residuals are computed with the signals from the 1801th minute to the 2000th minute. Each of the three residuals is designed to reject a fault affecting one of the three actuators. The residuals rejecting the faults of actuator 1, 2, and 3 are, respectively, plotted in the top, middle, and bottom pictures of Figure 2. The first residual is clearly smaller than the two others, indicating that the hypothesis of a fault affecting the first actuator is the most likely one, which corresponds to the actually simulated fault.

fig2
Figure 2: Fault isolation residuals. The time unit is the minute. Top: residual rejecting the fault of actuator 1. Middle: residual rejecting the fault of actuator 2. Bottom: residual rejecting the fault of actuator 3.

6. Conclusion

Despite unknown nonlinear distortions of sensors, the information provided by such sensors is still useful for fault diagnosis, even when there is no redundant sensors, if the distortions are strictly monotonous. The monotonousness is a weak assumption since nonmonotonous distortion would make the sensor information useless. The main idea of the method presented in this paper is about how to use the information provided by such sensors. Because of the unknown nature of nonlinear distortion, neither the absolute value of the measured physical variable nor its sign can be determined from the sensor signal. The strict monotonousness of the nonlinear distortion is not helpful in this aspect. However, for any two different time instants, the relative sign of the measured variable is preserved by the monotonous nonlinear distortion. By using the information residing in the relative sign of sensor signals, the method for actuator fault diagnosis presented in this paper is conceptually robust to sensor distortions, as illustrated by the numerical example presented in this paper.

Notes

1For notation simplicity, the sampling period is assumed to be 1 here.

References

  1. M. Basseville and I. Nikiforov, Detection of Abrupt Changes: Theory and Applications, Information and System Sciences Series, Prentice-Hall, Englewood Cliffs, NJ, USA, 1993.
  2. J. J. Gertler, Fault Detection and Diagnosis in Engineering Systems, Marcel Dekker, New York, NY, USA, 1998.
  3. J. Chen and R. J. Patton, Robust Model-Bsed Fault Diagnosis for Dynamic Systems, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1999.
  4. M.-A. Massoumnia, G. C. Verghese, and A. S. Willsky, “Failure detection and identification,” IEEE Transactions on Automatic Control, vol. 34, no. 3, pp. 316–321, 1989. View at Publisher · View at Google Scholar · View at MathSciNet
  5. C. De Persis and A. Isidori, “A geometric approach to nonlinear fault detection and isolation,” IEEE Transactions on Automatic Control, vol. 46, no. 6, pp. 853–865, 2001. View at Publisher · View at Google Scholar · View at MathSciNet
  6. P. W. Tse, W.-X. Yang, and H. Y. Tam, “Machine fault diagnosis through an effective exact wavelet analysis,” Journal of Sound and Vibration, vol. 277, no. 4-5, pp. 1005–1024, 2004. View at Publisher · View at Google Scholar
  7. H. Ma, C. Han, X. Kong, G. Wang, J. Xu, and X. Zhu, “A chaos-GFRF based fault diagnosis method,” in Proceedings of IEEE Conference on Cybernetics and Intelligent Systems (ICCIs '04), vol. 2, pp. 1314–1317, Singapore, December 2004. View at Publisher · View at Google Scholar
  8. K. Kumamaru, K. Inoue, and T. Iwamura, “Improved approach to KDI-based fault detection for non-linear black-box systems,” in Proceedings of the SICE Annual Conference, vol. 1, pp. 927–932, Sapporo, Japan, August 2004.
  9. M. E. H. Benbouzid, H. Nejjari, R. Beguenane, and M. Vieira, “Induction motor asymmetrical faults detection using advanced signal processing techniques,” IEEE Transactions on Energy Conversion, vol. 14, no. 2, pp. 147–152, 1999. View at Publisher · View at Google Scholar
  10. Q. Zhang, “A method for actuator fault diagnosis with robustness to sensor distortion,” in Proceedings of the 6th IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes (SAFEPROCESS '06), Beijing, China, August-September 2006.
  11. Q. Zhang, A. Iouditski, and L. Ljung, “Identification of Wiener system with monotonous nonlinearity,” in Proceedings of IFAC Symposium on System Identification, Newcastle, Australia, March 2006.
  12. R. B. Kearfott and V. Kreinovich, Eds., Applications of Interval Computations, R. B. Kearfott and V. Kreinovich, Eds., Kluwer Academic Publishers, Boston, Mass, USA, 1996.
  13. D. Den Hertog, Interior point approach to linear, quadratic, and convex programming: algorithms and complexity, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1994.
  14. T. Glad and L. Ljung, Control Theory: Multivariable and Nonlinear Methods, Taylor & Francis, New York, NY, USA, 2000.