Structured controller design for closed‐loop
 D
 ‐stability in convex/non‐convex regions: Mixed integer‐linear programming approach

Correspondence , Farshad Merrikh-Bayat, Department of Electrical and Computer Engineering, University of Zanjan, Zanjan, Iran. Email: f.bayat@znu.ac.ir Abstract Here it is assumed that the characteristic function of a linear control system, which is in the form of a polynomial and its coefficients are desired affine functions of unknown parameters of controller, is given. It is also assumed that the transfer function of controller has a desired order and structure. Hence, some of the coefficients of the characteristic polynomial may not depend on some or any of the controller parameters. The main problem under consideration is to calculate the parameters of the controller such that all roots of the characteristic equation lie (if possible) inside the desired D-stability region, which can be any convex/non-convex connected/disconnected subset of complex plane. This problem has important applications in control theory. For example, non-convex D-stability regions appear when designing a controller for a fractional-order system is aimed. It is shown that this problem, which is still open even in dealing with general convex regions, is exactly equivalent to a set of linear algebraic equalities/inequalities in mixed real-integer variables, which can be solved efficiently by using the available software. Application of the proposed method for optimal controller design is also studied and three numerical examples are presented.


INTRODUCTION
Stability is the crucial requirement of any control system. In practice, however, stability is not sufficient for good behaviour of the control system and it is desired that the closed-loop poles lie in a subset of complex plane which has some nice properties. For example, one may want to put a lower bound on both the exponential decay rate and the damping ratio of the closed-loop response. This requirement is fairly fulfilled by locating the roots of characteristic equation in a suitable region of complex plane. By definition, a polynomial is D-stable if all its roots are in the region D of the complex plane. A more rigorous mathematical definition is as follows. Let Γ be a closed Jordan curve and let int Γ and ext Γ denote its interior and exterior, respectively. Let D denote int Γ. A characteristic polynomial of degree n is D-stable if and only if the image of Γ under characteristic equation encircles the origin n times in a counter-clockwise manner [1,Ch. 3].

FIGURE 1
A non-convex D-stability region with desired properties for locating the closed-loop poles of a discrete-time system (FOS); see for example [7][8][9][10]. Because of the difficulty of taking into account the intrinsic non-convexity of the stability region in dealing with FOS, the few available studies on the D-stability of FOS are limited to convex regions, which lead to conservative results [11]. Non-convex regions also appear in studying the stability of discrete-time linear control systems. Assume that we want the settling time, rise time and damping ratio of a discrete-time system be equal to t s , t r and , respectively, or less. For this purpose we can first determine the suitable region for locating the closed-loop poles in s-plane and then use the mapping z = e sT , where T is the sampling period, to determine the corresponding region in z-plane [12,Ch. 7]. Using the formulas ℜ(s) = n ≥ 4.6∕t s , n ≥ 1.8∕t r [12, p. 213], and | arg(s)| ≥ − cos −1 and assuming t r = 1s, t s = 5s, = 0.6 and T = 0.1s the D-stability region in z-plane is obtained as shown in Figure 1. Note that in order to limit the consumption of control, the constraint ℜ(s) ≥ −6 is also considered in the s-plane, where ℜ(s) stands for the real part of s. Another application of D-stability for controller design is represented in the following. Consider a standard errorfeedback system where the transfer functions of controller and plant are equal to C (s) and P (s), respectively. Our aim is to design a controller which stabilises the closed-loop system and makes the closed-loop transfer function as close as possible to the desired model M (s). Ideally, such a controller satisfies the equation C (s)P (s) 1 + C (s)P (s) = M (s), (1) or equivalently, But, in general, for two main reasons the C (s) obtained from (2) is not acceptable. The first reason is that it does not have the desired structure while favourite controllers like proportionalintegral-derivative (PID) and lead-lag are low order with a specific structure. Recall that in the context of LTI systems, a controller is called full order and unstructured if its degree is equal to the degree of the corresponding plant to be controlled and any entry of its state-space matrices can be adjusted freely. The other reason is that it does not guarantee the stability of closedloop system. A possible remedy is to calculate C (s) from the following optimisation problem: where k are the samples taken from the frequency range of interest (totally L samples). Clearly, the controller obtained from (3) stabilises the closed-loop system and minimises the difference between the frequency response of the closed-loop system and the desired model. The optimisation problem in (3) equivalently can be written as follows: s.t. C (s) stabilises the feedback system, (4b) where ℜ and ℑ stand for real and imaginary part, respectively. When C (s) is linear in its unknown parameters (e.g. it has the structure of a PID), the inequality constraints in (4c) and (4d) are also linear in these unknown parameters. In this case (4a)-(4e) constitute a linear program which can be solved very effectively by the available software. But, unfortunately, the crucial stability constraint (4b) makes the problem much more complicated and harder to solve. In practice, for reasons like preventing the saturation of actuators, the D-stability constraint is considered instead of (4b). In Example 2 of Section 3 it is shown how the above problem can be formulated as a mixed integer-linear program (MILP) and then solved. Many of the methods developed in recent years for designing low-order structured controllers are based on solving a constrained optimisation problem by employing a gradient-based iterative algorithm. These methods usually begin from an initial stabilising controller and iteratively update it until a (local) optimal solution is found. For example, in [13] a weighted sum of the integral of the absolute errors caused by input and output disturbances is minimised subject to constrains on the peaks of the sensitivity and complementary sensitivity functions using the proposed exact gradients. A very similar idea is studied in [14]. In [15] parameters of a filtered PID are calculated by minimising the integral error (IE) performance index (in response to load disturbance) subject to constraints on the probability of violation of the H ∞ -norm of sensitivity and complementary sensitivity functions. Again the algorithm is gradient based and starts from an initial solution which stabilises the feedback system. In some other references like [16] and [17] the information of sub-gradients is used to solve a constrained H ∞ -norm minimisation problem. These methods are also of iterative nature and begin the search from a stabilising controller. The important point is that the solution obtained at the end by the abovementioned iterative algorithms is, in general, a local optimum whose quality drastically depends on the initial solution used to begin the search. In order to be confident that the structured controller obtained at end has some nice properties, the search should begin from a suitable D-stabilising controller.
Here the standard unity-feedback LTI system is considered and it is assumed that the characteristic equation of the closedloop system with a structured controller in the loop is in the form of ∑ n k=0 a k s k = 0 where some of a k 's are fixed real numbers and the others are affine functions of the parameters of controller. Using this assumption, an algorithm for calculating the parameters of controller such that all roots of the characteristic equation lie inside a pre-specified convex/non-convex D-stability region is proposed. This problem is still open even when the region is convex. A rather similar problem is studied in [18] where it is aimed to find a minimum order controller such that all poles of the MIMO feedback system are placed inside a pre-specified convex region. That study is based on the properties of polynomial matrices and has the advantage of being directly applicable to MIMO systems. However, the stability region under consideration is necessarily convex and the original non-convex problem is approximated by a convex one.
The method proposed for structured D-stabilising controller design finds the parameters of controller by solving an MILP. In recent years MILPs have found many applications in different fields of engineering such as train speed trajectory optimisation [19], optimisation of unit commitment and economic dispatch [20], optimal configuration planning for energy hubs [21], etc. Few applications of MILP in the field of control theory are also reported in the literature. In these applications first often a kind of linear (affine) mathematical model is obtained for the dynamical system under consideration by introducing some auxiliary integer (decision) variables. This approach is especially common in dealing with hybrid systems which involve both continuous and discrete variables. At the next step often a model predictive controller (MPC) is designed for this mathematical model of system by solving an MILP. For example, in [22] first a piecewise affine (PWA) system is modelled as a set of linear equalities and inequalities by introducing some binary and continuous-time auxiliary variables. Then it is shown that applying an explicit MPC (eMPC) to the resulting model yields a multi-parametric MILP problem which can be solved using the available software. As another example, in [23] first a microgrid management system is modelled by a set of mixed integer-linear equality/inequality constraints. Then these constraints are incorporated in an MPC for further compensation of uncertainties through the feedback mechanism. In this manner, at each time sample the optimal control law is obtained by solving an MILP. Paper [24] is also a classical reference on modelling complex systems by means of logical and integer decision variables. Although MILP is a very powerful tool for modelling complex systems and problems, it is known to be NP-complete, which means that the resulting program may be extremely difficult to solve in practice. Hence, special care should be taken during formulation of the problem to make sure that the resulting MILP is fairly relax and solvable by using the available commercial software.
The remaining is organised as follows. Main results, including the proposed MILP for designing a D-stabilising structured controller are presented in Section 2. Three numerical examples are studied in Section 3 and finally Section 4 concludes the paper.

Structured D-stabilising controller design
In this paper the characteristic equation of the feedback system is considered as f (s, k) = 0 where s = x + iy is the Laplace variable and k is the vector containing the unknown parameters of the (structured) controller under consideration. Clearly, for SISO systems f is always an affine function of k. For example, assuming the transfer functions of plant and controller as P (s) = 1∕(s + 1) and C (s) = (as + b)∕(s + c ), respectively, the characteristic equation can be written as f (s, k) = The aim of the following discussions is to develop a method for calculating a vector k, which contains the unknown parameters of the structured controller C (s), such that all roots of f (s, k) = 0 lie inside a desired convex/non-convex region in s-plane. Obviously, depending on the shape and size of the convex/non-convex region under consideration and the special structure considered for controller, this problem can have infinitely many or even no solutions. Hence, in general, the feasible k is not unique.
The first step of the proposed method for calculating a Dstabilising structured controller is meshing the region under consideration. For example, Figure 2 shows a non-convex region and a possible meshing where the nodes of mesh are identified by black spheres. As mentioned earlier, such a non-convex region which comprises points from the right half plane (RHP) appears in studying the stability of FOS. In this paper it is assumed that the nodes are located at the intersection points of some horizontal and vertical lines, and each cell is in the shape of a rectangle. As can be observed in Figure 2, only the cells whose all four nodes belong to the region under consideration are taken into account (these cells are identified by grey painting in Figure 2).
The above discussion turns out that when the border of region consists of slant lines, which is the case in Figure 2, the FIGURE 2 A non-convex region and a possible approximate meshing. No nodes should lie on the real axis cells only can approximately fill inside the region. However, the error caused by this approximation is not a big problem since one can make the difference between the approximate and the original region arbitrary small by making use of smaller cells. As another point, it will be shown later in this section that the width and length of the cells can be different. Moreover, in practice the meshing should be such that no nodes lie on the negative real axis. The reason is that the proposed algorithm puts the roots of the characteristic equation only inside the cells and if the negative real axis does not pass through some of the cells (i.e. the edges and vertices of some of the rectangles lie exactly on the negative real axis), the characteristic equation cannot have any negative real roots. This situation is impossible in dealing with many problems.
In the following, necessary and sufficient conditions to have a root inside a cell will be developed. For the sake of simplicity, first consider the parameter-free equation f (s) = u(x, y) + iv(x, y) = 0 where u = ℜ( f (s)) and v = ℑ( f (s)) are the real and imaginary parts of f (s), respectively. Obviously, the roots of this equation are the intersection points of the curves u(x, y) = 0 and v(x, y) = 0 in x − y plane. Suppose that u(x, y) = 0 and v(x, y) = 0 have an intersection point (of order one) at s 0 = x 0 + iy 0 which is located inside a rectangular cell as shown in each of the sub-figures of Figures 3 and 4.
Considering the fact that u is a continuous two-variable function of x and y, it must be positive for the points located on one side of the curve u(x, y) = 0 and negative for the points located on its other side. A same discussion goes on v(x, y). Depending on the relative positions of the curves u(x, y) = 0 and v(x, y) = 0 with respect to each other, and the signs of u and v on the two sides of these curves, totally exactly eight scenarios are possible as depicted in Figures 3 and 4. On the other hand, when f is an analytic function, which is the case in our problem since f is a polynomial, the Cauchy-Riemann necessary conditions for analyticity of the characteristic polynomial must also be sat- hold at this point. But the situations shown in Figure 4 violate the Cauchy-Riemann necessary conditions for analyticity at s = s 0 , which means that only the situations shown in Figure 3 can happen in practice (e.g. it is observed in Figure 4a that u x (x 0 , y 0 ) < 0 and v y (x 0 , y 0 ) > 0 which concludes that u x (x 0 , y 0 ) ≠ v y (x 0 , y 0 ) and the Cauchy-Riemann necessary conditions are violated).
In the problem under consideration f is the characteristic function of the closed-loop system which depends on k. Hence, according to the above discussion, equation f (s, k) = u + iv = 0 has a root (of order one) inside the rectangular cell under consideration if and only if the signs of u and v satisfy exactly one of the cases shown in Figure 3. This statement is equivalent to saying that exactly one of the following four logical sentences must be true: where ∧ stands for the logical AND, and the points s 11 , s 12 , s 21 and s 22 are considered as shown in Figure 5. The main problem with the logical sentences in (6) is that one cannot directly solve them to find a feasible k. The other problem is that the number of the roots of f (s, k) = 0 in the region of interest (e.g. as the one shown in Figure 2) is equal to the degree of f (s, k), while logical trueness of exactly one of the sentences in (6) guarantees placing only one of the roots of f (s, k) = 0 in one of the cells of the mesh of interest. In the following, by introducing suitable decision variables we transform the sentences in (6) to a set of algebraic relations in k whose solution puts each root of f (s, k) = 0 in one of the cells of the mesh under consideration. For the sake of simplicity in notation, f (s, k) is shown as f (s) in the following discussions.
In order to get rid of the logical ANDs in (6) and arrive at equivalent algebraic relations, assign five binary decision variables to each cell (denoted as 1 , 2 , 3 , 4 and 1 ) and two binary decision variables to each node of the mesh (denoted as R 11 , R 12 , R 21 , R 22 , I 11 , I 12 , I 21 and I 22 ) as shown in Figure 5. These binary decision variables are defined as the following: Using the above definition for the binary decision variables, for each cell of the mesh write five relations as follows: where is a sufficiently small positive real constant (e.g. = 0.01). Considering the fact that for any binary variable X we have X = 1 − X , inequalities in (15)- (18) can be written as follows: Now, for any jkth node (i.e. the node located at the j th row and the kth column) of the mesh write four mixed integerbinary inequalities as follows where m and M are sufficiently large negative and sufficiently large positive real constants (e.g. m = −10 4 and M = 10 4 ). Note that in the right-hand side of inequalities (23)-(26) R jk and I jk are unknown binary decision variables, and on the left-hand side of these inequalities the vector k (which is not shown for the sake of simplicity) contains unknown real variables. Hence, relations (23)-(26) are actually mixed real-binary inequalities.
To sum up, for any cell of the mesh write one equation in binary decision variables of the cell as equation (14), and eight inequalities in binary decision variables of the cell and the nodes of that cell as given in relations (19)- (22), and for any node of the mesh write four mixed real-binary inequalities as given in relations (23)- (26). At the end, add the following equality in binary variables to the set of existing relations where N is the total number of cells in the mesh under consideration, n is the degree of the characteristic function, and i is the binary variable assigned to the ith cell (recall that each cell has one and four variables which are related as given in Equation 14). Considering the fact that i 's in Equation (27) are binary variables, it is concluded from this equation that exactly n out of the totally N , i variables will be equal to unity and the others will be equal to zero. In the following it is shown that when i = 1 then exactly one root of the characteristic equation sits inside the ith cell. Note that before solving the corresponding MILP it is not known for sure which i 's will be equated to unity and which cells are selected by the algorithm to put the roots of characteristic equation inside them. Of course, it may also happen that the problem has no solutions. For example, in Figure 5 assuming 1 = 1 and considering the fact that i 's (i = 1, … , 4) are binary variables, it is concluded from Equation (14) that only one of the i 's in this equation can be equal to one and the others must be equal to zero. The i which is equal to one is determined only after solving the program. For now, without any loss of generality, assume that 1 = 1 and 2 = 3 = 4 = 0. Substitution of these values in inequalities (15) The binary variables cannot be determined uniquely from inequalities (29)-(31), but inequality (28) For m → −∞, M → ∞ and ↓ 0 one can easily verify that inequalities in (32)-(35) are exactly equivalent with the conditions represented in sentence (6a), as is expected. It is concluded from the above discussion that one solution for equation 1 = 1 + 2 + 3 + 4 = 1 is 1 = 1 and 2 = 3 = 4 = 0 which implies trueness of sentence (6a). Similarly, it can be verified that another solution for equation 1 = 1 + 2 + 3 + 4 = 1 is 2 = 1 and 1 = 3 = 4 = 0 which implies trueness of sentence (6c), and so on. Hence, by setting 1 = 1 the algorithm makes one of the i 's equal to one and the others equal to zero, which makes exactly one of the sentences in (6a)-(6d) logically true and others logically false (note to the definition of i 's in the conditional Equations (10)- (13)). Recall that each of the cases represented in sentences (6a)-(6d) address one of the possible cases for f (s, k) = 0 to have a root inside the cell identified with the binary variable i . Figure 6 shows a sample mesh and the corresponding binary decision variables assigned to its nodes and cells. As can be observed in this figure, a mesh with 4 cells and 9 nodes requires 20 binary variables for cells and 18 binary variables for nodes.
Hence, the binary variables required to solve a problem which consists of many cells can be very large. In practice, we can mesh the D-stability region under consideration using a few cells and solve the problem (the number of cells must be greater than or equal to the degree of the closed-loop characteristic equation). Then, if the algorithm failed to find a solution, use more number of smaller cells to mesh the region and solve the problem again. Repeat this procedure until a solution is found (provided that at least one solution exists).

Summary of the proposed algorithm
where a 1 = −1.5, a 2 = −0.5, b 1 = −0.15, and b 2 = 1. With these values for parameters the process has one RHP pole and one non-minimum phase (NMP) zero. Hence, unavoidable initial undershoot and limitations on performance are expected. The controller under consideration is an FOPD with transfer function: which is located in series with process in a standard errorfeedback loop. In this example the time constant of the derivative filter, , is considered equal to 0.05. Recall that, considering the characteristic equation as f (s) = 1 + C (s)P (s) = 0, this closed-loop system is bounded-input bounded-output (BIBO) stable if and only if all roots off (w) := 1 + C (s)P (s)| s =w = 0 lie inside the non-convex sector defined by [7,8]: In this example two types of problems can be solved using the proposed method: the feasibility problem and the gain optimisation problem. The aim of the feasibility problem is to calculate k p and k d in (37) such that all roots off (w) = 0 lie inside the non-convex D-stability region shown in Figure 7, where the border of region is identified by solid lines. Note that this D-stability region is a subset of ℂ 0.5 . The gain optimisation problem is defined as maximising |k p | subject to the constraint that all roots off (w) = 0 lie inside the D-stability region under consideration. The solution obtained for the gain optimisation problem minimises the tracking and disturbance rejection errors at steady state. For this purpose we can solve both min k p (39) and min −k p (42) where D is the set of all points which belong to the D-stability region. After solving the above two problems, we consider the k p with the larger amplitude and the corresponding k d as the final solutions. First assume = 0.5. By employing the step-by-step algorithm of Section 2.2 with the meshing shown in Figure 7, the gains of the optimal controller are obtained as k p = 20.6848 and k d = 6.7588. In this case the corresponding MILP consists of 2202 variables which takes about 83 s to solve using MAT-LAB R2014b running on an i5 2.20-GHz Intel core processor with 8GB of RAM. Locations of the closed-loop poles with this optimal FOPD controller are identified by × in Figure 7. Interesting point is that the algorithm has placed closed-loop poles in the corners of the non-convex D-stability region which shows the importance of having access to an algorithm which can deal with non-convex regions. Solving the feasibility problem yields k p = 0.7448 and k d = 5.3092. The corresponding closed-loop poles are identified by + in Figure 7 in this case. The MILP of the feasibility problem is solved much faster than the optimal controller design problem, approximately in 4.26s.
For the sake of simplicity, in order to simulate the response of the closed-loop system to a step command we assume that both the controller and plant are initially at rest, which means that all initial conditions are equal to zero. Under this assumption, the well-known definitions for fractional derivative like Riemann-Liouville, Grünwald-Letnikov and Caputo lead to the same transfer functions [25,Ch. 2]. Then, the time-domain step response of the closed-loop system is obtained by calculating the numerical inverse Laplace transform of the FO closed-loop transfer function multiplied by 1∕s using the method developed in [26]. The Matlab function for numerical evaluation of inverse The response of the closed-loop system to the unit step command is shown in Figure 8 for four different cases. As is expected and also observed in this figure, the controller obtained by solving the gain maximisation problem for = 1∕2 minimises the steady-state error compared to the controller obtained by solving the feasibility problem. This improvement in steady-state error is obtained at the cost of larger overshoots in the step response (compare the dash-dotted plot with the dotted one in Figure 8). The reason is the optimal controller locates the complex-conjugate closed-loop poles closer to the border of instability region compared to the feasible controller (compare the location of the poles identified as × with those identified as + in Figure 7). No considerable overshoot is observed in the simulations of Figure 8 corresponding to = 1∕3. The reason is that in this case the sector of stability in w-plane is | arg(w)| > ∕2 = ∕6, and consequently, the closed-loop poles are fairly faraway from the border of instability. Recall that some initial undershoot in the response is unavoidable because of the NMP zero of process.
Example 2. Consider again the FO process of Example 1 as defined in (36) and assume a 1 = −1.5, a 2 = 0.5, b 1 = 0.5, b 2 = 10, and = 0.5. With these values for parameters the process has two unstable poles at 1 and 0.25. Our aim here is to design an FOPI controller with transfer function such that firstly all of the closed-loop poles lie in the D-stability region of Figure 9, and secondly, frequency response of the closed-loop system is as close as possible to the frequency response of the desired model defined as This problem can be formulated as a constrained optimisation problem as discussed in Section 1, and then solved using the proposed method. Note that in the problem under consideration, constraints (4c) and (4d) yield: and where C k , P k and M k stand for C ( j k ), P ( j k ) and M ( j k ), respectively, and L is the number of frequency samples. Hence, the main difference between the mathematical formulation of this problem and the problem studied in the previous example is that here we must add extra constraints in the form of (4e), (47) and (48) to the problem and also consider the cost function as given in (4a). In this problem the frequency range of interest is considered as 0 < ≤ 100 rad/s and 150 logarithmically spaced frequency samples are taken from it (i.e. 0 < k ≤ 100 rad∕s, k = 1, … , 150). Figure 10 shows the response of the feedback system (with the optimal controller in the loop) and the corresponding desired model to the unit step command. In this figure the controllers are designed assuming two different models for FIGURE 10 Response of the closed-loop system and desired model to unit step command, corresponding to Example 2

FIGURE 11
Effect of the size of cells on the value of cost function and the optimal values of k p and k i , corresponding to Example 2 the desired closed-loop system: one with = 0.5 and n = 10 rad/s, which yields k p = 0.5501 and k i = 0.5055, and the other with = 0.5 and n = 100 rad/s, which yields k p = 2.3234 and k i = 8.3878. As it is observed in this figure, step response of the closed-loop system changes effectively according to the changes made in the transfer function of model. Poles of the closed-loop system with the resulting two different controllers are shown in Figure 9. In this figure + and × corresponds to the models with n = 100 rad/s and n = 10 rad/s, respectively. Figure 11 shows the effect of the size of cells on the value of cost function as defined in (4a) and the optimal values obtained for k p and k i . In order to calculate the data points in this figure the non-convex D-stability region is considered similar to Figure 9 and it is assumed that the nodes have a same vertical and horizontal distance from each other. Then the side of square cells is varied between 0.5 and 2 with steps of 0.1, and in each case the values of cost function, k p and k i are calculated. As is observed, the plot of cost function versus the length of sides is fairly flat which means that the algorithm is not so sensitive to the size of cells.
Example 3. Consider a two-input two-output process with the following state-space matrices: which has four poles at 1.2065, 0.2493, −3.0238 and −2.0320. It can easily be verified that this system is controllable, which means that all of the closed-loop poles can be placed at the predetermined desired locations by applying a full static state feedback. But, here it is assumed that the state-feedback matrix is structured as follows: and the aim is to calculate k 11 , k 13 , k 22 ∈ ℝ (if possible) such that all poles of the closed-loop system lie inside the rectangle defined as {s ∈ ℂ : −8 < |ℜ(s)| < −3, |ℑ(s)| < 5.5}, in s-plane. Note that by applying the state-feedback matrix (51) we use a limited information of states to generate controls (state number four is not used at all and the second control is calculated only based on the information of the second state). Figure 12 shows the D-stability region as defined in (52), the nodes of the mesh under consideration and the closedloop poles after solving the problem. The state-feedback gains are obtained as k 11 = 21.5078, k 13 = 1.5573 × 10 4 and k 22 = 1.4562 × 10 4 . Clearly, by modifying the state-space matrices, or structure of the state-feedback matrix, or the definition of D-stability region we may arrive at a problem with no solutions.

CONCLUSION
A method for designing a structurally constrained controller to achieve closed-loop D-stability in a desired convex/non-convex region is developed. It was especially shown that this problem is exactly equivalent to finding a solution to a set of linear algebraic equalities/inequalities in mixed integer/linear variables. Although this problem is NP-complete, the proposed formulation leads to fairly relax equalities/inequalities, which can be solved effectively by a trivial software like Matlab. In fact, problems with thousands of mixed integer/linear variables are solved successfully and the results are reported. Three numerical example were also designed and solved by using the proposed method. In two of these examples the Dstability region is non-convex and the designed optimal controller locates the poles of closed-loop system in the corners of non-convex region. This observation justifies the importance of non-convex corners and having access to an algorithm that can deal with them. In the third example it was shown that the proposed method can be used to calculate the static state-feedback gains of a two-input two-output process such that all poles of the closed-loop system lie inside the pre-specified D-stability region. This problem is also challenging since the state-feedback matrix is assumed to be structured.
One limitation of the proposed method is that all poles of the resulting closed-loop system will be necessarily distinct. In fact, taking into account the possibility of occurring repeated poles drastically increases the complexity of the problem due to the new decision variables required for its formulation. However, it is not a big loss of generality because the distinct poles can appear arbitrarily close to each other and mimic the behaviour of repeated poles.