1 Introduction

Feedback controls are most popular in industry [1]. A great deal of effort has been made in designing optimal feedback control gains for various applications [2]. In the time domain, for example, the overshoot, rise time or peak time, settling time and the tracking error are often used to characterize the performance of the closed loop system. It is well-known that the overshoot and peak time are conflicting objectives, meaning that when the overshoot goes down, the peak time goes up, and vice versa. It is thus quite natural to consider the multi-objective feedback control design to minimize the overshoot, peak time, and tracking error at the same time. The multi-objective optimization algorithm can find a set of optimal solutions that represent the best compromises among these conflicting goals. This paper presents the simple cell mapping (SCM) method for multi-objective optimal design of feedback controls for linear dynamical systems with or without time delay.

There have been many multi-objective optimal control studies in the literature [3]. Rani et al. [4] developed the global ranking genetic algorithm to design a proportional–integral–derivative (PID) controller for a rotary inverted pendulum. Three conflicting functions are considered in the design including the settling time, overshoot and mean square steady-state error. Another multi-objective genetic algorithm called the genetic artificial immune system algorithm was investigated by Khoie et al. [5]. The four fitness functions are considered: the rise time, overshoot, settling time, and integral square error. A non-dominated sorting genetic algorithm was used by Kumar and Nair [6] to tune the control gains to minimize the rise time, overshoot and settling time. Similar multi-objective optimal control designs are studied with different algorithms for multi-objective optimization problems (MOPs) including the ant colony algorithm by Ibtissem et al. [7], an immune algorithm by Kim [2] and the imperialist competitive algorithm by Esmaeil and Caro [8]. Other designs are in the frequency domain where the gain and phase margins are used to characterize the system stability and robustness, and the crossover frequency is used to assess the system response speed [9]. Liu and Daley [10] have proposed an optimal PID control for a rotary hydraulic system. The control parameters are tuned such that the crossover frequency, gain margin, phase margin, and steady state error are within the targeted range.

The solution of MOPs does not consist of a single point in the design space, but rather forms a set, called Pareto set named after the Italian economist Vilfredo Pareto (1848–1923) [11]. The corresponding objective function values are called Pareto front. Many numerical methods for MOPs have been studied. There exist, for instance, scalarization methods that transform the MOP into a scalar optimization problem (SOP). By choosing a clever sequence of SOPs, a finite size approximation of the entire Pareto set can be obtained in certain cases [12–15]. See the book by Eichfelder [16] for an overview. These methods are advantageous in particular for uni-modal (e.g., convex) objectives due to their fast convergence when classical mathematical programming techniques are used to solve the SOPs. However, they may run into trouble for multi-modal objectives due to the locality of the approach. It is possible that these methods get stuck in local optimal solutions that are not globally optimal. Since the Pareto set forms under some mild regularity conditions locally a manifold, the continuation methods which perform searches along the Pareto set are very efficient if one solution is at hand and if the Pareto set is connected [17–20].

Evolutionary algorithms are most widely used for MOPs [21]. The underlying idea in evolutionary computation is to steer (or evolve) an entire set of solutions (population) toward the set of interest during the search. Evolutionary algorithms specialized for multi-objective problems have been shown their strength in many applications. Due to the global and stochastic nature of evolutionary algorithms, the Pareto set can be approximated quite well in most cases, although there is always an uncertainty left regarding whether the global Pareto set has indeed been found. There also exist attempts to compute nearly optimal solutions with evolutionary algorithms [22]. However, this approach suffers the drawback that only a subset of the nearly optimal solutions can be stored. In particular, a significant fraction of the set of interest may be ignored by the archiving strategy in certain cases [22].

Another approach to approximate the Pareto set is to use the set oriented methods with subdivision techniques [23–25].The advantage of the set oriented methods is that they generate an approximation of the global Pareto set in one single run of the algorithm. Further, they are applicable to a wide range of optimization problems and are characterized by a great robustness. Hence, these methods are interesting alternatives against ‘classical’ mathematical programming techniques in particular for the thorough investigation of low or moderate dimensional MOPs. The cell mapping method in this study is the predecessor of the set oriented methods, and was proposed by Hsu [26] for global analysis of nonlinear dynamical systems. In the cell mapping method for MOPs, the dynamical systems are derived from multi-objective optimization search algorithms. The cell mapping method can obtain a finite size approximation \(\mathcal{A }\) of the Pareto set \(\mathcal{P }\) such that the distance between \(\mathcal{P }\) and \(\mathcal{A }\) is less or equal to a given threshold value in the Hausdorff sense. That is, the approximation quality of \(\mathcal{A }\) can in principle be determined. A first implementation of the multi-objective cell mapping technique is used that already yields promising results.

Two cell mapping methods have been extensively studied, namely, the SCM and the generalized cell mapping (GCM) to study the global dynamics of nonlinear systems [26, 27]. The cell mapping methods have been applied to optimal control problems of deterministic and stochastic dynamic systems [28–30]. Other interesting applications of the cell mapping methods include optimal space craft momentum unloading by Flashner and Burns [31], single and multiple manipulators of robots by Zhu and Leu [32], optimum trajectory planning in robotic systems by Wang and Lever [33], and tracking control of the read-write head of computer hard disks by Yen [34]. Sun and his group studied the fixed final state optimal control problems with the SCM method [35, 36] , and applied the cell mapping methods to the optimal control of deterministic systems described by Bellman’s principle of optimality [37]. Crespo and Sun further applied the generalized cell mapping based on the short-time Gaussian approximation to stochastic optimal control problems [30, 38]. They also studied semi-active optimal control of populations of competing species in a closed environment with the cell mapping method [39]. This paper for the first time applies the SCM method to study multi-objective optimal control problems.

In this paper, we consider two time domain design problems of feedback controls for linear systems with or without time delay to demonstrate the SCM method for MOPs. In Sect. 2, we present the linear control systems. In Sect. 3, we describe the MOP formulation and the SCM method applied to the MOP. Then, the MOP of the feedback control design is presented with the goal to simultaneously minimize the peak time, overshoot and integrated absolute tracking error. In Sect. 4, we present numerical examples of the two linear control systems, one of which has a control time delay, and consider the design in the feedback gain space and with the help of \(\mathbf{Q}\) and \(\mathbf{R}\) matrices in the linear quadratic regulator (LQR) optimal control formulation.

2 Linear control system

Consider a general linear control system

$$\begin{aligned} {\dot{\mathbf{x}}}(t)&= \mathbf{Ax}(t)+\mathbf{Bu}(t), \\ \mathbf{y}(t)&= \mathbf{Cx}(t)\nonumber \end{aligned}$$
(1)

where \(\mathbf{x}\in \mathbf{R}^{n},\; \mathbf{u}\in \mathbf{R}^{m}\) and \(\mathbf{y}\in \mathbf{R}^{p}\). Consider a full state feedback control \(\mathbf{u}=-\mathbf{K}\mathbf{x}(t)\). The control gain \(\mathbf{K}\) can be designed in a number of ways. When \((\mathbf{A,B})\) is a controllable pair, we can design the control by the method of pole placement. For single-input–single-output (SISO) systems, the control can also be designed in the time domain by considering the time domain specifications or in the frequency domain by considering the stability margins [40]. One of the popular optimal control design methods is the LQR for both SISO and multi-input–multi-output (MIMO) systems. The control is found to minimize a quadratic cost function,

$$\begin{aligned} J=\int \limits _{0}^{\infty }(\mathbf{x}^{T}\mathbf{Qx}+ \mathbf{u}^{T}\mathbf{Ru})dt \end{aligned}$$
(2)

subject to the constraint of the state equation as well as initial and terminal conditions. The matrices \(\mathbf{Q}\) and \(\mathbf{R}\) determine the relative weights of the response and the control effort, and ultimately determine the control gain. The sub-optimal control gain is given by \(\mathbf{K}=\mathbf{R}^{-1}\mathbf{B}^{T}\mathbf{P}\) where \(\mathbf{P}\) satisfies the algebraic Riccati equation [41]. Note that for tracking problems, the LQR control can be formulated in terms of the tracking error.

2.1 Time-delayed linear control system

Consider now the system with a control time delay \(\tau \).

$$\begin{aligned} {\dot{\mathbf{x}}}(t)&= \mathbf{Ax}(t)+\mathbf{Bu}(t-\tau ), \\ \mathbf{y}(t)&= \mathbf{Cx}(t).\nonumber \end{aligned}$$
(3)

We define an extended state vector as

$$\begin{aligned} \mathbf{z}(t)\!=\!\left[ \mathbf{x}(t),\;\mathbf{u}(t\!-\!\tau _{N}),\;\mathbf{u}(t\!-\!\tau _{N-1}),\cdots ,\mathbf{u}(t\!-\!\tau _{1})\right] ^{T},\quad \end{aligned}$$
(4)

where \(\mathbf{z}\in \mathbf{R}^{n+mN}\), \(0=\tau _{0}<\tau _{1}<\cdots <\tau _{N-1}<\tau _{N}=\tau \), and \(N\) is the number of mesh grids to discretize the interval \([0,\tau ] \). After introducing an interpolation scheme of \(\mathbf{u}(t-\tau _{i})\) over the mesh grid \(\tau _{i}\) (\(i\in [0,N] \)), we can obtain an equation for \(\mathbf{z}(t)\) without an explicit time delay of the control to replace Eq. (3) as

$$\begin{aligned} {\dot{\mathbf{z}}}(t)&= {\bar{\mathbf{A}}z}(t)+{\bar{\mathbf{B}}u} (t),\\ \mathbf{y}(t)&= {\bar{\mathbf{C}}z}(t).\nonumber \end{aligned}$$
(5)

We should point out that Eq. (5) is now in the same form as Eq. (1). The control design can now be done in the extended state space \(\mathbf{R}^{n+mN}\) using the same methods as discussed above.

Note that when the system contains a retarded element as a function of the delayed response \(\mathbf{x}(t-\tau _{s})\), the extended state vector can be defined in the same way by including the delayed response \(\mathbf{x} (t-\tau _{j})\) where \(\tau _{j}\) is the mesh grid on \([ 0,\tau _{s}] \), \(j\in [ 0,N_{s}] \), and \(N_{s}\) is the number of mesh grids to discretize the interval \(\left[ 0,\tau _{s}\right] \). This method is known as continuous time approximation. It can handle multiple independent time delays of the system and the control, and accurately predict a wide range of infinite poles of the time-delayed system. More details of the method can be found in [42, 43].

3 Multi-objective optimization

A multi-objective optimization problem (MOP) can be expressed as follows:

$$\begin{aligned} \min _{\mathbf{k}\in Q}\{\mathbf{F}(\mathbf{k})\}, \end{aligned}$$
(6)

where \(\mathbf{F}\) is the map that consists of the objective functions \(f_{i}:Q\rightarrow \mathbf{R}^{1}\) under consideration.

$$\begin{aligned} \mathbf{F}:Q\rightarrow \mathbf{R}^{k},\; \mathbf{F}(\mathbf{k})=[f_{1} (\mathbf{k}),\ldots ,f_{k}(\mathbf{k})]. \end{aligned}$$
(7)

\(\mathbf{k}\in Q\) is a \(q\)-dimensional vector of design parameters. The domain \(Q\subset \mathbf{R}^{q}\) can in general be expressed by inequality and equality constraints:

$$\begin{aligned} Q&= \{\mathbf{k}\in \mathbf{R}^{q}\ |\ g_{i}(\mathbf{k})\le 0,\quad i=1,\ldots ,l,\\&\quad \mathrm{and}\quad h_{j}(\mathbf{k})=0,\quad j=1,\ldots ,m\}.\nonumber \end{aligned}$$
(8)

Next, we define optimal solutions of a given MOP by using the concept of dominance [11].

Definition 1

  1. (a)

    Let \(\mathbf{V},\mathbf{W}\in \mathbf{R}^{k}\). The vector \(\mathbf{V}\) is said to be less than \(\mathbf{W}\) (in short: \(\mathbf{V}<_{p}\mathbf{W}\)), if \(V_{i}<W_{i}\) for all \(i\in \{1,\ldots ,k\}\). The relation \(\le _{p}\) is defined analogously.

  2. (b)

    A vector \(\mathbf{v}\in Q\) is called dominated by a vector \(\mathbf{w}\in Q\) (\(\mathbf{w}\prec \mathbf{v}\)) with respect to the MOP (6) if \(\mathbf{F}(\mathbf{w})\le _{p}\mathbf{F}(\mathbf{v})\) and \(\mathbf{F}(\mathbf{w})\ne \mathbf{F}(\mathbf{v})\), else \(\mathbf{v}\) is called non-dominated by \(\mathbf{w}\).

If a vector \(\mathbf{w}\) dominates a vector \(\mathbf{v}\), then \(\mathbf{w}\) can be considered to be a ‘better’ solution of the MOP. The definition of optimality or the ‘best’ solution of the MOP is now straightforward.

Definition 2

  1. (a)

    A point \(\mathbf{w}\in Q\) is called Pareto optimal or a Pareto point of the MOP (6) if there is no \(\mathbf{v}\in Q\) which dominates \(\mathbf{w}\).

  2. (b)

    The set of all Pareto optimal solutions is called the Pareto set denoted as

    $$\begin{aligned} \mathcal{P }:=\{\mathbf{w}\in Q:\mathbf{w}\quad {\mathrm{is}}\; {\mathrm{a}}\; {\mathrm{Pareto}}\; {\mathrm{point}}\; {\mathrm{of}}\; {\mathrm{the}}\; {\mathrm{MOP}}\; (6) \}.\nonumber \\ \end{aligned}$$
    (9)
  3. (c)

    The image \(\mathbf{F}(\mathcal{P })\) of \(\mathcal{P }\) is called the Pareto front.

Pareto set and Pareto front typically form (\(k-1\))-dimensional manifolds under certain mild assumptions on the MOP. See [18] for a thorough discussion.

3.1 SCM method

The cell mapping methods describe system dynamics with cell-to-cell mappings by discretizing both the state space and time. It starts with a point-to-point mapping as a finite difference approximation of the governing differential equation of the system as

$$\begin{aligned} \mathbf{x}(k)=\mathbf{G}(\mathbf{x}(k-1)), \end{aligned}$$
(10)

where \(\mathbf{x}\left( k\right) \in \mathbf{R}^{n}\) is the state vector at the \(k\)th mapping step. Then, the SCM proposes to discretize the state space into a set of small cells, and represents the dynamics of an entire cell denoted as \(Z\) by the dynamics of its center. The center of \(Z\) is mapped according to the point-to-point mapping. The cell that contains the image point is called the image cell of \(Z\). The cell-to-cell mapping is denoted by \(C\),

$$\begin{aligned} Z(k)=C(Z(k-1)). \end{aligned}$$
(11)

We should note that the exact image of the center of \(Z\) is approximated by the center of its image cell according to the SCM. This approximation can cause significant errors in the long term solution of dynamical systems [28, 29, 33]. Nevertheless, the SCM offers an effective approach to investigate global response properties of the system.

The cell mapping with a finite number of cells in the computational domain will eventually lead to closed groups of cells of the period same as the number of cells in the group. The periodic cells represent invariant sets, which can be periodic motion and stable attractors of the system. The rest of the cells form the domains of attraction of the invariant sets. For more discussions on the cell mapping methods, their properties and computational algorithms, the reader is referred to the book by Hsu [26].

To apply cell mapping techniques to compute MOP solutions, we need to define dynamical systems derived from the search algorithms. As will be shown below, the invariant set of the dynamical system represents the Pareto set of the MOP.

Assume that the design parameter \(\mathbf{k}\) is updated by a map \({\varvec{\Phi }}\) that generates descent directions of all the objective functions \(f_{i}\) at a given point \(\mathbf{k}\in Q\) [13, 44, 45]. Let \({\varvec{{\nu }}}\in \mathbf{R}^{q}\) be a descent direction, along which a dominating solution can be found. That is, for \(\mathbf{k}\) and \({\varvec{\nu }}\) there exists a \(t_{0}\in \mathbf{R}_{+}\) such that \(\mathbf{k} \,+\, t_{0} {\varvec{\nu }}\prec \mathbf{k}\). Assume that the map is defined by the following equation,

$$\begin{aligned} \mathbf{k}(n)&= \mathbf{k}(n-1)+{\varvec{\Phi }}(\mathbf k (n-1)),\\ \mathbf{k}(0)&= \mathbf{k}_{0}\in Q.\nonumber \end{aligned}$$
(12)

The descent map \({\varvec{\Phi }}\) with a suitable step size can be used to search for the invariant set in \(Q\). In the following, we present an example of the descent map by Fliege and Svaiter [13]. Define an auxiliary function,

$$\begin{aligned}&\mathbf{g}:\mathbf{R}^{q}\rightarrow \mathbf{R}^{k}\nonumber \\&\mathbf{g}({\varvec{\nu }}) = \max _{i=1,\ldots ,k}\left[ \mathbf{J} (\mathbf{k}){\varvec{\nu }}\right] _{i}, \end{aligned}$$
(13)

where \(\mathbf{J}(\mathbf{k})\) denotes the Jacobian of \(\mathbf{F}\) at \(\mathbf{k}\). \(\mathbf{g}\) is convex and positive homogeneous. With the help of \(\mathbf{g}\), the following convex optimization problem defines a descent map,

$$\begin{aligned} {\varvec{\Phi }}_{F}(\mathbf{k})=\min _{{\varvec{\nu }}\in \mathbf{R}^{k}}\left[ \mathbf{g}({\varvec{\nu }})+\frac{1}{2}\Vert {\varvec{\nu }} \Vert _{2}^{2}\right] . \end{aligned}$$
(14)

Note that the Jacobian is calculated over the cell partition in the framework of cell mapping. Therefore, it is a finite difference approximation of the Jacobian. We should also point out that other search algorithms will lead to different maps. For a comprehensive survey of various search algorithms, see the book by Liu et al. [46].

3.2 Multi-objective optimal design of feedback controls

Now, consider the multi-objective optimization approach to design feedback control gains of linear dynamical systems with or without time delay. In general, different objective (fitness) functions can be considered for optimization in time domain or frequency domain. For instance, peak time and overshoot are considered in time domain, while phase and gain margins are considered in frequency domain [1, 10, 47]. Here, we consider the multi-objective optimization design in time domain.

In particular, we consider the following MOP to design the control gain \(\mathbf{k}\) or the weighting matrices \(\mathbf{Q}\) and \(\mathbf{R}\) in the performance index of the LQR control (2),

$$\begin{aligned} \min _{\mathbf{k}\in Q}\left\{ t_{p},M_{p},e_{IAE}\right\} , \end{aligned}$$
(15)

where \(Q\) represents the set of admissible control gains \(\mathbf{k}\) or the set of the weighting matrices \(\mathbf{Q}\) and \(\mathbf{R}\), \(M_{p}\) stands for the overshoot of the response to a step reference input, \(t_{p}\) is the corresponding peak time and \(e_{IAE}\) is the integrated absolute tracking error

$$\begin{aligned} e_{IAE}=\int \limits _{0}^{T_{ss}}\left| r(\hat{t})-x(\hat{t})\right| d\hat{t}. \end{aligned}$$
(16)

where \(r(t)\) is a reference input and \(T_{ss}\) is the time when the response is close to be in the steady state. When the control magnitude is finite, the two factors, i.e. peak time and overshoot, are conflicting, while the integrated tracking error plays a compromising role.

4 Numerical examples

4.1 First order system plus time delay

Consider a first order plus time delay (FOPTD) system,

$$\begin{aligned} X(s)=\frac{K}{T_s+1}e^{-\tau s}U(s), \end{aligned}$$
(17)

where \(K\) and \(T\) are constants and \(\tau \) is the time delay. Define a tracking error as \(e=x-r\) where \(r\) is a reference input. We augment the plant with an additional state \(x_{2}\) such that \(\overset{.}{x}_{2}=x-r\), i.e. \(x_{2}= {\displaystyle \int _{0}^{t}} e(t)dt\). We set \(x_{1}=x\) and \(\mathbf{x}=[x_{1},x_{2}]^{T}\). The augmented state equations read

$$\begin{aligned} \left\{ \begin{array} [c]{c} \overset{.}{x_{1}}\\ \overset{.}{x_{2}} \end{array} \right\} \!&= \!\left[ \begin{array} [c]{ll} -\frac{1}{T} &{}\quad 0\\ 1 &{}\quad 0 \end{array} \right] \left\{ \begin{array} [c]{l} x_{1}\\ x_{2} \end{array} \right\} \!+\!\left[ \begin{array} [c]{l} \frac{K}{T}\\ 0 \end{array} \right] u(t\!-\!\tau ) \!+\!\left[ \begin{array} [c]{l} 0\\ -1 \end{array} \right] r(t). \nonumber \\ \end{aligned}$$
(18)

Note that a comparison study between the classical PI controller and the fractional order PI for this system is done by Luoa and Chen [48]. Here, we first convert the system (18) to the extended state equation (5) without time delay, and then digitize it. We then apply the discrete time LQR formation of the system (5) to design the full extended state feedback control as

$$\begin{aligned} u(k)&= -\,\mathbf{Kz}(k)=-\mathbf{k}_{x}\mathbf{x}(k)-k_{1}u(k-N) \nonumber \\&-\,k_{2} u(k-N+1),\ldots ,k_{N}u(k-1), \end{aligned}$$
(19)

where

$$\begin{aligned} \mathbf{z}(k)&= \left[ \mathbf{x}(k),u(k\!-\!N),u(k\!-\!N\!+\!1),\ldots ,u(k\!-\!1)\right] ^{T},\nonumber \\ \mathbf{x}(k)&= \mathbf{x}(k\Delta t),\\ u(k-i)&= u(k\Delta t\!-\!\tau _{i}),\quad i\in \left[ 0,N\right] ,\nonumber \\ \mathbf{K}&= [\mathbf{k}_{x},k_{1},\ldots ,k_{N}],\nonumber \end{aligned}$$
(20)

\(\mathbf{k}_{x}\) is \(1\times 2\) and represents the PI control gain for the original time-delayed system. \(\Delta t\) denotes the sample time. The performance index of the LQR optimal control reads

$$\begin{aligned} J=\frac{1}{2} {\displaystyle \sum \limits _{k=0}^{\infty }} \mathbf{z}^{T}(k)\mathbf{Qz}(k)+u(k)Ru(k). \end{aligned}$$
(21)

We have chosen the matrix \(\mathbf{Q}=diag(Q_{1},Q_{2},0,\ldots ,0)\) in this example. In general, the weighting on the delayed control can be non-zero.

In the numerical example reported next, we have chosen \(\tau =0.5\) s, \(K=1\) and \(T=1\). The sample time of the digital control is \(\Delta t=0.05\) s. The discretization number of the time delay is \(N=10\). We take \(\mathbf{k} =[Q_{1},\,Q_{2},R]\) as the design parameters of the MOP. The design space for the parameters is chosen as follows,

$$\begin{aligned}&Q=\left\{ [Q_{1},Q_{2},R]\in [1,3]\times [200,500]\right. \nonumber \\&\qquad \left. \times [0.011,0.0134]\subset \mathbf{R}^{3}\right\} . \end{aligned}$$
(22)

In the SCM method, we select the number of partitions of the design space \(Q\) as \(\mathbf{N}=[20,30,10]\) without further subdivisions. The total CPU time for computing all the solutions is \(27\) min on a laptop PC.

Figure 1 shows the Pareto set of the MOP solution consisting of 1,198 cells and Fig. 2 shows the corresponding Pareto front. In the LQR control formulation, different parameters \([Q_{1},\,Q_{2},R]\) with the same proportionality will lead to the same optimal feedback gain \(\mathbf{K}\) and therefore the same Pareto front points. This explains why the Pareto front in Fig. 2 contains so much overlap as compared to the Pareto front of the next example. The peak time of the Pareto front in Fig. 2 is uniformly spaced at the interval equal to the sample time of the digital control. To yield a more dense distribution of the peak time, one can adopt continuous time control design instead of digital control. It is interesting to examine the closed-loop poles of the system on the Pareto set. Since the LQR design is done in the extended state space, \(N\) gains \([k_{1},\ldots ,k_{N}]\) form a filter to regulate the history of the control. The gain vector \(\mathbf{k}_{x}\) describes the original system dynamics. We plot the closed-loop poles of the original system (18) on the Pareto set in Fig. 3. The dot on the right represents the integrator pole and the line segment on the left contains the poles of the first order system. Figure 4 shows the closed-loop step response and the corresponding control with the gains that lead to the smallest \(e_{IAE}\) in the Pareto set. The control performance is quite satisfactory.

Fig. 1
figure 1

The Pareto set of \([Q_{1},Q_{2},R]\) for the multi-objective LQR optimal control of the first order system with time delay. The color code indicates the level of the other design variable. Red denotes the highest value, and dark blue denotes the smallest value. (Color figure online)

Full size image
Fig. 2
figure 2

The Pareto front of \([t_{p},M_{p},e_{IAE}]\) for the multi-objective LQR optimal control of the first order system with time delay corresponding to the Pareto set in Fig. 1. The color code indicates the level of the other objective function. Red denotes the highest value, and dark blue denotes the smallest value. (Color figure online)

Full size image
Fig. 3
figure 3

The first two dominant closed-loop poles for the multi-objective LQR optimal control of the first order system with time delay corresponding to the Pareto set in Fig. 1. The color code indicates the level of the first objective function \(t_{p}\). Red denotes the highest value, and dark blue denotes the smallest value. (Color figure online)

Full size image
Fig. 4
figure 4

An example of the closed-loop step response and control of the first order system with time delay. The PID gains are designed with the multi-objective LQR optimal control formulation. The optimal design parameters are \([Q_{1},Q_{2},R]=[3.10,495.0,0.012]\). The optimal gains for the original state \(\mathbf{x}\) are \(\mathbf{k}_{x}=[53.62,110.93]\). The performance indices are \([t_{p},M_{p},e_{IAE}]=[1.05,0.1839,0.4407]\).

Full size image

4.2 Second order linear oscillator

Consider a second order oscillator subject to a PID control.

$$\begin{aligned} \ddot{x}+2\zeta \omega _{n}\dot{x}+\omega _{n}^{2}x=\omega _{n}^{2}u(t), \end{aligned}$$
(23)

where \(\omega _{n}=5\), \(\zeta =0.01\),

$$\begin{aligned} u(t)=k_{p}\left[ r(t)-x(t)\right] +k_{i}\int \limits _{0}^{t}\left[ r(\hat{t}\,)-x(\hat{t}\,)\right] d\hat{t}-k_{d}\dot{x}(t),\nonumber \\ \end{aligned}$$
(24)

\(r(t)\) is a step input, \(k_{p}\), \(k_{i}\) and \(k_{d}\) are the PID control gains. We consider the MOP in Sect. 3.2 with the control gains \(\mathbf{k}=\left[ k_{p},k_{i},k_{d}\right] ^{T}\) as design parameters. The design space for the parameters is chosen as follows,

$$\begin{aligned} Q=\left\{ \mathbf{k}\in [10,50]\times [1,30] \times [1,2]\subset \mathbf{R}^{3}\right\} . \end{aligned}$$
(25)

Initially, we select the number of divisions in the three control gain intervals as \(\mathbf{N}=[30,18,8]\). The integrated absolute tracking error \(e_{IAE}\) is calculated over time with \(T_{ss}=20\) seconds. After the first run of the SCM program, we refine the Pareto set with \(7\times 7\times 7\) subdivisions. The total CPU time for this example is \(20\) min. The closed-loop response of the system for each design trial is computed with the help of closed form solutions.

Figure 5 shows the Pareto set of the MOP solution and Fig. 6 shows the corresponding Pareto front. 5,961 cells are in the Pareto set. In this case, the Pareto front exhibits a fine structure, which has not been seen before, and cannot be readily obtained with random search algorithms such as the genetic algorithm. Figure 7 shows the closed-loop poles of the system on the Pareto set. The general location of the cluster of the poles is along the \(\pm 45^{\circ }\) lines to the left of the imaginary axis, which is consistent with the well-known feedback control intuitions [40]. Figure 8 shows the closed-loop step response and the control with a gain that leads to the smallest \(e_{IAE}\) in the Pareto set. The performance functions are indicated in the figure caption.

Fig. 5
figure 5

The Pareto set of \([k_{p},k_{i},k_{d}]\) for the multi-objective optimal control of the second order system. The color code indicates the level of the other design variable. Red denotes the highest value, and dark blue denotes the smallest value

Full size image
Fig. 6
figure 6

The Pareto front of \([t_{p},M_{p},e_{IAE}]\) for the multi-objective optimal control of the second order system corresponding to the Pareto set in Fig. 5. The color code indicates the level of the other objective function. Red denotes the highest value, and dark blue denotes the smallest value

Full size image
Fig. 7
figure 7

The closed-loop poles for the multi-objective optimal control of the second order system corresponding to the Pareto set in Fig. 5. The color code indicates the level of the second objective function \(M_{p}\). Red denotes the highest value, and dark blue denotes the smallest value.

Full size image
Fig. 8
figure 8

An example of the closed-loop step response and control of the second order system with \([k_{p},k_{i},k_{d}]=[40.0,2.8796,1.9792]\). The PID gains are designed with the multi-objective optimal control. \([t_{p},\,M_{p},\,e_{IAE}]=[0.1555,0.0,0.2774]\)

Full size image

5 Concluding remarks

We have reviewed the MOP formulation and the SCM method applied to the MOP. We have then studied the multi-objective optimal time domain design of feedback controls for linear systems with or without time delay with the help of the SCM method. It should be pointed out that we have only considered simple feedback control examples for the purpose of demonstrating the proposed SCM–MOP design method. There are different structures of controls that render perfect tracking. We have considered two different sets of design parameters for the feedback control: a LQR based approach with the weighting matrices as design parameters, and a direct optimization with feedback gains as design parameters. Both approaches prove to be quite effective. The Pareto set and Pareto front consisting of the peak time, overshoot and integrated absolute tracking error are obtained for examples of two linear control systems, one of which has a control time delay. It is interesting to note that for the second order linear system, we have found a fine structure of the Pareto front, which has been very difficult to obtain using stochastic search algorithms such as the genetic algorithm. This study suggests that the SCM method is a powerful method that can provide global and fine-structured solutions of MOPs for complex dynamical systems.