Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2010 (2010), Article ID 105143, 20 pages
Research Article

The Design of QFT Robust Compensators with Magnitude and Phase Specifications

1Departamento de Lenguajes y Computación, Escuela Politécnica Superior, Universidad de Almeria, 04120 Almeria, Spain
2Departamento de Informática y Sistemas, Facultad de Informática, Universidad de Murcia, 30100 Murcia, Spain

Received 27 April 2010; Revised 29 October 2010; Accepted 27 December 2010

Academic Editor: Jerzy Warminski

Copyright © 2010 José Carlos Moreno et al. 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.


The frequency response is an important tool for practical and efficient design of control systems. Control techniques based on frequency response are of special interest to dealing with important subjects such as the bandwidth and the cost of feedback. Furthermore, these techniques are easily adapted to deal with the uncertainty of the process to control. Quantitative feedback theory (QFT) is an engineering design technique of uncertain feedback systems that uses frequency domain specifications. This paper analyzes the phase specifications problem in frequency domain using QFT. This type of specification is not commonly taken into account due to the fundamental limitations of the linear control given by Bode's integral. An algorithm is proposed aimed at achieving prespecified closed-loop transfer function phase and magnitude variations, taking into account the plant uncertainty. A two-degrees-of-freedom feedback control structure is used and a new type of boundary is defined to satisfy these objectives. As the control effort heavily depends on a good estimation of these boundaries, the proposed algorithm allows avoiding overdesign.

1. Introduction

Feedback around a plant is mandatory because of uncertainty, including uncertain disturbances acting on . The main task of feedback in control must be to achieve a desired input-output relation within desired tolerances despite plant uncertainty. However, not all control techniques deal with the uncertainty in an explicit form; in particular in the literature the control techniques dealing with the uncertainty are called robust control techniques [1]. Another fundamental subject in feedback control is the quantitative sense. The plant uncertainty and the desired tolerances of the input-output relation must be able to be quantitatively formulated.

Feedback is fundamental in order to achieve a particular behavior of a system but unfortunately there exist limitations. The greater the bandwidth of a system, the faster the tracking of reference input, but the sensitivity to sensor noise is augmented. Also, the frequency range over which one must be reasonably confident of the uncertainty of the system is increased. This is the cost of the feedback [2], and it is important to be able to deal with it in a quantitative form.

The benefits and cost of feedback are strongly frequency dependent, so design in this domain has the added advantage of direct, powerful control of the system's feedback properties and costs [2]. Also, the trade-off between compensator complexity and bandwidth economy is highly visible in this domain.

Quantitative Feedback Theory is a frequency domain robust control technique that can be considered as a framework to design robust controllers, where the system uncertainty is typically of parametric nature, commonly given in form of templates (see [3] for an introduction to QFT). QFT uses a two-degrees-of-freedom control scheme (see Figure 1), where the uncertain system is represented by a transfer function belonging to a family of plants , and where and are the compensator and the precompensator, respectively, which must be designed in order to achieve performance and stability robust specifications.

Figure 1: The typical 2DoF control scheme used in QFT.

In QFT, the specifications on input-output relations are given in the frequency domain, in terms of admissible bounds on the frequency response of the closed-loop transfer functions between the different inputs and outputs. For a number of design frequencies, these specifications are combined with the system uncertainty description (templates) in order to obtain constraints, usually referred to as boundaries in the QFT literature. For each design frequency , boundaries are given as curves in the Nichols plane, delimiting allowable regions for (or equivalently for , where is the nominal plant). In a second step, nominal specifications are used to shape the precompensator .

This paper focuses on the analysis of the problems associated with the simultaneous consideration of magnitude and phase specifications for the closed loop transfer function from reference input to system output (which can be of interest in high-precision tracking problems as for instance coordinated movement in robotics, CNC, flight simulators, and in dealing with nonlinear elements, as actuator saturation, in the control loop, as will be pointed out later).

Few works about phase specifications and their applications can be found in the QFT literature [49]. The way in which these papers address the problem differs from the approach used in this paper, which mainly consists of shaping and to achieve some nominal phase and magnitude specifications using a new set of boundaries. Another important point considered is the computation of multiple-valued boundaries. For instance, the algorithm proposed in [5] to compute the phase tracking boundaries did not exploit the fact that boundaries can be multiplevalued.

The consideration of multiple-valued boundaries may have an important practical relevance, as the control effort is directly related to them. This fact was pointed out in [10] and considered in the subsequent works, but general solutions to this problem have not been found. The computation of multiple-valued tracking boundaries has been analyzed in [11] and extended in [12] to include phase tracking boundaries to guarantee certain phase variations on closed loop transfer function from reference input to system output. In this work, a new type of boundary (the nominal phase tracking boundary) is proposed to cope with the nominal phase specifications.

As has been pointed out, phase specifications appear in some practical applications, for example in coordinated motion, and have been studied in QFT in several works including [5, 8]. In robotics, the robot motion control problems can be separated into two categories: positioning and contouring. In contouring problems, the robot tool tip is commanded to follow a specific path. Here, the spatial contour tracking accuracy of the robot is of paramount concern, since it directly influences the quality of the final product. In other contour control systems like the computerized numerical control (CNC) machines, contouring accuracy is also crucial. All of these cases can be handled using phase specifications. This type of specifications are also important for flight simulators and in general for any high-accuracy tracking system. Feed-forward controllers like zero phase error tracking controllers (ZPETC) [13] and cross-coupled controllers (CCC) have been developed to effectively reduce tracking error and contouring error, respectively. They have been used in some cases with robust controllers, usually or QFT, to take the plant uncertainty into account [1416]. Other papers in this topic are based on a feedforward control for set-point tracking and QFT [1719]. It is important to note that in [14, 1619] the precompensator in Figure 1 is always fixed to 1, and the uncertainty is taken into account to design compensator only for the magnitude specifications. Other authors deal with the phase in an implicit form using the magnitude of the tracking error [6, 7, 9] however; the study of phase specifications is important by itself. For example, in [20], phase specifications are used in the design of feedback systems with actuator saturation in order to guarantee the stability of a nonlinear control system.

The design of control systems with phase specifications (and uncertainty) is hard to handle, basically because it is an optimization problem with nonconvex constraints, a problem without a closed solution. In this paper, the QFT framework is used to deal with the problem. The phase constraints are transformed into an additional set of boundaries, and the QFT methodology is used to solve the controller optimization/design problem.

The paper is organized as follows. After some preliminaries in Section 2, the subsequent sections show different approaches for solving the phase specification problem. Section 3 considers a first solution for a particular case, when the plant is an integrator with uncertain gain. In Section 4, a general solution is investigated. An example is developed in Section 5, and, finally, conclusions are outlined in Section 6.

2. Preliminaries

The problem that will be considered in this work is the design of a control system (Figure 1) to satisfy tracking specifications, considering a nominal value of the closed loop transfer function from reference input to system output , (the complex variable is omitted in general for simplicity in the notation), and allowed deviations. is any element of a family of LTI (linear time invariant) plants , and is the nominal plant. For the nominal value and allowed variations, the specifications are given for both magnitude and phase.

can be designed to meet variations over both magnitude and phase of the closed loop transfer function . The role of the precompensator is to fix the nominal value of , , but due to the fact that phase and magnitude of are related, by Bode's integral assuming minimum phase systems, or analogous constraints for unstable and nonminimum phase systems, phase and magnitude can not be independently manipulated in design.

Usually, and are designed without taking into account the phase specifications in the design process, then a satisfactory design can be obtained meeting variations over the nominal magnitude of , but this is not the general case.

In the following, some key definitions for the design method are introduced.

Definition 2.1. Let with the open loop transfer function and the complementary sensitivity function of system in Figure 1, then the following sets , with the set of plants defining the border of the typical template in QFT, for each frequency and are defined.

Definition 2.2. is a crossing set of functions if and for some frequency . The rank of frequencies in which there are crossings is noted by , where . Equivalently a set of functions is noncrossing if with .

Definition 2.3. For a crossing set of functions , is defined as , the maximum difference between the magnitud of transfer functions belonging to for all frequencies in .

In this paper tracking specifications given by

In [7] the tracking specifications are given in form of relative error on sensitivity function as shown in.

In [4, 9] a similar type of specification is used. In Figure 2, the three most well-known approaches to deal with typical tracking specifications for lowmedium frequencies are shown: specifications in classical QFT, specifications from [7], and specifications given by (2.1) and (2.2).

Figure 2: Comparison between the three most well-known forms of dealing with the tracking specifications in QFT.

3. A Solution for a Particular Case

In order to satisfy the set of closed-loop specifications in (2.1) and (2.2), the first step is, given design specifications in (2.4) and (2.5) over a frequencies set (Figures 3 and 4), to find in order to satisfy them. Once has been obtained, fulfilling the corresponding set of bounds obtained from the specifications and the templates, has to be designed to achieve nominal specifications in (2.1) and (2.2) (Figures 5 and 6). In [12] an algorithm to compute the boundaries over the nominal open-loop transfer function (equivalently over ) is presented. This algorithm is based on the construction of a 3-D surface, where the boundaries are simply contour lines, and here is used to compute the tracking magnitude boundaries for specifications in (2.4), and the tracking phase boundaries for specifications in (2.5). After boundaries computation, the compensator is shaped to fulfill them, so that the specifications in (2.4) and (2.5) are satisfied, and, finally, the precompensator is designed.

Figure 3: Specifications on the variation of magnitude of in  (2.4).
Figure 4: Specifications on the variation of phase of in  (2.5).
Figure 5: Specifications on the magnitude of .
Figure 6: Specifications on the phase of .

The main problem in the design of is that, in our case, two objectives (nominal phase and magnitude specifications in (2.1) and (2.2)) have to be met with only one degree of freedom, . A first approach to cope with this problem is summarized in Algorithm 3.1.

Algorithm 3.1   3.1. We have the following steps.Step 1   1. Shape the compensator fulfilling boundaries for specifications in (2.4) and (2.5), computed with the algorithm in [12]. Step 2   2. For each frequency compute such that taking Step 3   3. Choose .

In Step 2, the point has been calculated for each , in such a way that the M-contour passing through the point is the maximum M-contour passing through the . If is a noncrossing set of functions then the computation of is very simple, due to the fact that maximum M-contour passes through the same point () for all , so that . In general, the algorithm guarantees the fulfillment of the specifications in (2.1) but not necessarily the satisfaction of specifications in (2.2) (in general, the maximum N-contour passing through each does not pass through the point ). Furthermore, in order to guarantee a proper selection of , the pole-zero excess of must be less than or equal to the pole-zero excess of . This last constraint can be overcome if the satisfaction of phase specifications is not demanded at high frequency (the normal situation in practical cases).

If the plant is a single integrator with uncertain gain, and a robust stability specification less or equal to 0 dB is used, then the maximum M and N contours pass through the same point of each . In this particular case, the above algorithm is a good solution for the problem. For example, if the uncertain plant and the compensator are given, respectively, by then a non-crossing set of open-loop transfer functions (Figure 7) and a crossing set (Figure 8) are obtained. So, there exists a value of in (3.3) such that the maximum in the set of open loop transfer functions is achieved for this value in all frequencies, but this is not the case when obtaining , due to the location of open loop transfer functions on the Nichols Plane (NP) and the form of the M-contour (an M-contour is the locus of points such that ) at this location. As shown in Figure 9, the member of set with maximum magnitude depends on the frequency, so that a crossing set as shown in Figure 8 is obtained.

Figure 7: Frequency response in Nichols plane for a noncrossing set of open-loop transfer functions computed from (3.3) and (3.4).
Figure 8: Bode diagrams for the functions belonging to the crossing set computed from (3.3) and (3.4).
Figure 9: Justification of crossings in set.

The specifications for this example are given by the following functions:

In this case, the tracking specifications given by (2.4) and (2.5), and a robust stability specification of 0.2 dB are used, in order to show that an adequate approximation to the tracking specifications in (2.1) and (2.2) can be achieved, using Algorithm 3.1.

The following set of design frequencies rps is considered. The tracking specifications, from (2.4) and (2.5), for functions in (3.5) are shown in Table 1.

Table 1: Magnitude and phase tracking specifications.

Using the algorithm in [12], the magnitude and phase tracking boundaries are computed (Figures 10 and 11). The nominal open loop transfer function (Figure 12) is obtained using the computer tools in [21]. The designed controller is given by (3.6), where a high-frequency pole is added in order to achieve a pole-zero excess of one. In Figure 13, the Bode diagrams for magnitude and phase of functions in the set , with given by (3.6) and given by (3.3), and specifications given by (3.5), are shown,

Figure 10: Magnitude tracking boundaries.
Figure 11: Phase tracking boundaries.
Figure 12: Nominal open-loop shaping.
Figure 13: transfer functions (-) and specifications (- -).

Then Algorithm 3.1 is applied to obtain that after some simplifications is given by

In Figure 14, magnitude and phase of the closed-loop transfer functions are shown approximately satisfying the specifications in (3.5). Function has been selected belonging to the set , but it can be observed in Figure 15 that this selection is incorrect, since is a crossing set of transfer functions. The approximation in this case has provided good results because is small. Obviously, there exists a relation between and the situation of nominal open loop transfer function in NP. So, when the is near point (,0 dB), is higher than when it is far. In this example, a robust stability specification of 0.2 dB has been used, this being related with the approximation error incurred when selecting belonging to the set .

Figure 14: transfer functions (-) and specifications (- -).
Figure 15: Set at low-medium frequencies for given by (3.3) and given by (3.6).

This example has shown that the computation of , in some cases, can be performed by defining an approximation error, given by a robust stability specification, and choosing as a function belonging to set . The approximation will be better as the plant is closer to an integrator.

4. A New Type of Boundary

The idea developed in this section is based on using and in Figure 1, in order to achieve a design satisfying nominal magnitude and phase specifications, given by (2.1) and (2.2). A new type of specification and a new type of boundary that must fulfill are considered instead of using only the precompensator for this purpose after is designed taking into account only the specifications in (2.4) and (2.5).

In the nominal open loop transfer function shaping, the allowed region of the NP has to be restricted to a zone such that maximum M- and N-contours passing through the cross the same point. This is possible only for a template with all its points lying at the same phase (an example is the plant in (3.3)). In general, a new specification has to be defined directly related with the approximation error. So, a new type of boundary is proposed in order to satisfy this new type of specification, which denominated the nominal phase shaping boundary. This boundary provides an allowed region , for nominal open loop transfer function in NP, given by the following expression with and being the nominal plant.

In order to compute this region, the algorithm in [12] is easily adapted. For each frequency , the template is shifted (denoted by over the NP so that, for each phase and each magnitude, the difference between the maximum N-contour passing through the and the N-contour passing through the point given by (4.2) (the maximum M-contour passes through this point) is computed. This generates a surface in a three-dimensional space, in such a way that the new boundary for frequency can be computed by taking a section of this surface corresponding to a constant value .

This algorithm can be formulated as the set of the following three steps.

Algorithm 4.1. We have the following steps.Step 1. Choose an array of phases and an array of magnitudes .Step 2. For each phase in and for each magnitude in :
if then
Step 3. Boundary = level curve of at the height of the specification . Where and represents the border.

The allowed region for shaping with this type of boundaries is the region above them if they are open, or the outer region if they are closed. Due to the shape of the N-contours in the NP, it can be assured that the fulfillment of this constraint is compatible with the fulfillment of constraints given by the magnitude and phase tracking boundaries in (2.1) and (2.2).

Furthermore, the shape of this new type of boundary can be characterised from the shape of the templates, which can be typified within two types and (see Figure 16).(i): the largest magnitude points are located at the greatest phase (right part of the template).(ii): the largest magnitude points are not located at the greatest phase.

Figure 16: Examples of the different types of templates.

The shape of the nominal phase shaping boundary will be:(i)closed (the allowed region is the outside of the boundary) for type templates. Due to the fact that at small magnitudes, the phase of the template points is the same that the N-contours passing through these points.(ii)open (the allowed region is above the boundary) for type templates. Obviously for large enough specifications the boundaries for this type of templates may be closed.

In addition, can be used in the algorithm as a parameter to obtain less conservative results (restrictiveless boundaries for greater values of . So, the design procedure to compute and in Figure 1 consists of the application of Algorithm 3.1, and updating Step 1 so that the new type of boundaries computed with Algorithm 4.1 are taken into account. It is important to note that Steps 2 and 3 in Algorithm 3.1 can be modified to use and instead of and , respectively. In this case, in (4.5) would have to be easily updated to adapt to this change. It is important to note that these results are applicable to minimum phase plants, nonminimum phase plants, and unstable plants.

5. An Illustrative Example

Consider the uncertain plant (taken from [22]) given by the array of frequencies , used for design purposes, is rps, the nominal plant is given by (5.1) with , and the tracking specifications in (2.1) and (2.2) are given by the following bounds

Firstly, magnitude and phase specifications in (2.4) and (2.5) are considered. Table 2 gives values of and for the design frequencies. A value of degrees is used as the nominal phase shaping specification for all frequency . Using Algorithm 4.1 to compute the nominal phase shaping boundaries and the algorithm in [12] in order to compute the tracking phase and magnitude boundaries, results shown in Figure 17 are obtained. Using computer tools [21] the nominal open loop transfer function is designed (Figure 17(d)). In this case, the compensator is given by where a high-frequency pole is added in order to achieve a pole-zero excess of one.

Table 2: Magnitude specifications for the tracking problem.
Figure 17: Phase (a): magnitude, (b): nominal phase shaping, (c): tracking boundaries and the nominal open loop transfer function fulfilling all of them (d).

Figure 18, containing the final result of the design stage using the proposed algorithm, shows how specifications are satisfied almost over the whole frequency axis. Specifications are not satisfied in interval , so a new boundary could be computed for  rps for example and should be redesigned introducing poles and zeros around this frequency, but in this case the redesign of the compensator is not necessary due to the fact that the specifications are not satisfied at very high frequency. The selected precompensator is given by

Figure 18: Variations over the magnitude and the phase of due to uncertainty in plant.

Suppose that only the magnitude tracking specifications from (2.1) are taken into account, as done traditionally in QFT, then the controller in (5.5) is computed to satisfy specifications in (2.4) for all frequency in set , fulfilling all the constrains as shown in Figure 19. After that, the precompensator, given by (5.6), is computed, using three notch filters to deal with the specifications at medium frequencies, to satisfy the specifications in (2.1). Figure 20 shows how the magnitude tracking specifications are fulfilled but not the phase specifications. So, the new set of specifications and corresponding boundaries must be taken into account if specifications on phase are required where denotes and denotes .

Figure 19: Magnitude tracking boundaries and the nominal open loop transfer function fulfilling all of them.
Figure 20: Variations over the magnitude and the phase of due to uncertainty in plant without considering phase specifications.

6. Conclusions

In this paper, the problem of dealing with phase specifications in frequency domain to achieve robust designs for uncertain systems has been studied, and an algorithm based on QFT has been proposed to solve it. The algorithm is based on the inclusion of a new types of boundary in the loop-shaping stage. The different approaches to deal with phase specifications are shown and compared. By means of the classical control structure used in QFT, a 2DoF controller is synthesized to solve the problem. A typical example in the QFT literature is used to compare the Bode diagram of a classical design without taking the phase specifications into account and the Bode diagram of a design resulting of the application of the new algorithm.


This paper has been supported by MCYT under grants DPI2007-66455-C02-01 and DPI2007-66718-C04-04.


  1. M. Sidi, Design of Control Systems: from classical to modern practical approaches, Krieger, Malabar, Fla, USA, 2001.
  2. M. I. Horowitz, Quantitative Feedback Design Theory (QFT), QFT Publications, Boulder, Colo, USA, 1993.
  3. I. Horowitz, “Quantitative feedback theory,” IEE Proceedings D, vol. 129, no. 6, pp. 215–226, 1982. View at Google Scholar · View at Scopus
  4. S. M. Mahdi Alavi, A. Khaki-Sedigh, B. Labibi, and M. J. Hayes, “Improved multivariable quantitative feedback design for tracking error specifications,” IET Control Theory and Applications, vol. 1, no. 4, pp. 1046–1053, 2007. View at Publisher · View at Google Scholar · View at Scopus
  5. F. N. Bailey and M. Kallel, “Loop gain-phase shaping techniques for robust multi-axis coordinated motion control,” in Proceedings of the American Control Conference, pp. 921–923, June 1992. View at Scopus
  6. E. Boje, “Pre-filter design for tracking error specifications in QFT,” International Journal of Robust and Nonlinear Control, vol. 13, no. 7, pp. 637–642, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  7. E. Eitelberg, “Quantitative feedback design for tracking error tolerance,” Automatica, vol. 36, no. 2, pp. 319–326, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  8. T. Holt and M. Y. Lee, “A control strategy for multi-axis robot contouring accuracy enhancement,” in Proceedings of 1st National Conference on Applied Mechanism and Robotics, vol. 3, Cincinnati, Ohio, USA, 1989. View at Scopus
  9. R. A. Perez, O. D. I. Nwokah, and D. F. Thompson, “Almost decoupling by quantitative feedback theory,” in Proceedings of the American Control Conference, pp. 2001–2006, Boston, Mass,USA, June 1991. View at Scopus
  10. F. N. Bailey and C. H. Hui, “Cacsd tools for loop gain-phase shaping design of siso robust controllers,” in Proceedings of IEEE Control System Society Workshop on Computer Aided Control System, pp. 151–167, Cincinnati, Ohio, USA, 1989. View at Scopus
  11. R. A. Perez, A. Baños, and F. J. Montoya, “An algorithm for computing qft multiple-valued performance bounds,” in Proceedings of 3th International Symposium on QFT and other Frequency Domain Methods and Applications, pp. 29–32, Stratchclyde, Scotland, 1997. View at Scopus
  12. J. C. Moreno, A. Baños, and M. Berenguel, “Improvements on the computation of boundaries in QFT,” International Journal of Robust and Nonlinear Control, vol. 16, no. 12, pp. 575–597, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  13. M. Tomizuka, “Zero phase error tracking algorithm for digital control,” Journal of Dynamic Systems, Measurement, and Control,, vol. 109, no. 4, pp. 65–68, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  14. Z. Guo, J. Cao, and K. Zhao, “Cross-coupled approach synchronization controller design method based on quantitative feedback theory for dual hydraulic motors driven flight simulators,” in Proceedings of the International Conference on Computer Design and Applications (IC-CDA '10), vol. 3, pp. 376–381, 2010. View at Scopus
  15. L. Wang, F. Jin, and Y. Sun, “Contour control for direct drive xy table,” in Proceedings of the IEEE International Conference on Mechatronics and Automation, pp. 4919–4923, Changchun, China, August 2009. View at Scopus
  16. S. S. Yeh and P. L. Hsu, “Analysis and design of the integrated controller for precise motion systems,” IEEE Transactions on Control Systems Technology, vol. 7, no. 6, pp. 706–717, 1999. View at Publisher · View at Google Scholar · View at Scopus
  17. Z. Guo, C. Wang, and K. Zhao, “Equal-status approach synchronization controller design method based on quantitative feedback theory for dual hydraulic driven flight simulators,” in Proceedings of the International Conference on Computer Design and Applications (ICCDA '10), vol. 3, pp. 69–74, 2010. View at Scopus
  18. Z. Guo, C. Wang, and K. Zhao, “Qft compound controller design method for improving workspace motion fidelity of hydraulic flight simulators,” in Proceedings of the International Conference on Information and Automation, June 2010. View at Scopus
  19. G. Zhifu, C. Jian, and Z. Keding, “Design and experiments of model-free compound controller of flight simulator,” Chinese Journal of Aeronautics, vol. 22, no. 6, pp. 644–648, 2009. View at Publisher · View at Google Scholar · View at Scopus
  20. J. C. Moreno, A. Baños, and M. Berenguel, “A synthesis theory for uncertain linear systems with saturation,” in Proceedings of the 4th IFAC Symposium on Robust Control Design, Milan, Italy, 2003. View at Scopus
  21. C. Borguesani, Y. Chait, and O. Yaniv, The Quantitative Feedback Theory Toolbox for MATLAB, The MathWorks, Natick, Mass, USA, 1995.
  22. I. M. Horowitz and M. Sidi, “Synthesis of feedback systems with large plant ignorance for prescribed time-domain tolerances,” International Journal of Control, vol. 16, no. 2, pp. 287–309, 1972. View at Google Scholar · View at Zentralblatt MATH · View at Scopus