Open AccessArticle

On Application of the Ray-Shooting Method for LQR via Static-Output-Feedback

Yossi Peretz

Department of Computer Sciences, Lev Academic Center, Jerusalem College of Technology, P.O. Box 16031, 93721 Jerusalem, Israel

Algorithms 2018, 11(1), 8; https://doi.org/10.3390/a11010008

Submission received: 29 October 2017 / Revised: 27 December 2017 / Accepted: 4 January 2018 / Published: 16 January 2018

(This article belongs to the Special Issue Algorithms for Hard Problems: Approximation and Parameterization)

Download Versions Notes

Abstract

In this article we suggest a randomized algorithm for the LQR (Linear Quadratic Regulator) optimal-control problem via static-output-feedback. The suggested algorithm is based on the recently introduced randomized optimization method called the Ray-Shooting Method that efficiently solves the global minimization problem of continuous functions over compact non-convex unconnected regions. The algorithm presented here is a randomized algorithm with a proof of convergence in probability. Its practical implementation has good performance in terms of the quality of controllers obtained and the percentage of success.

Keywords:

optimal control; state-space models; randomized algorithms; continuous-time systems

1. Introduction

Let a continuous-time system be given by

\{\begin{matrix} \dot{x} (t) = A x (t) + B u (t) \\ y (t) = C x (t) \end{matrix}

(1)

where

A \in R^{p \times p}, B \in R^{p \times q}, C \in R^{r \times p}

and

x, u, y

are the state, the input and the measurement, respectively. We assume that

(A, B)

and

(A^{T}, C^{T})

are controllable (see [1] for a new reduction from the case where

(A, B)

and

(A^{T}, C^{T})

are only stabilizable). Let

J (x_{0}, u) : = \int_{0}^{\infty} (x {(t)}^{T} Q x (t) + u {(t)}^{T} R u (t)) d t,

(2)

denote the cost functional, where

Q > 0

and

R \geq 0

. Assuming that

x_{0} = x (0)

is given, the LQR problem is to attenuate the disturbance

x_{0}

using minimal control cost, i.e., to design a regulation input

u (t)

that minimizes

J (x_{0}, u)

, subject to the constraints given by (1). Let

u = - K y

be the static-output-feedback (SOF) with the closed-loop matrix

A_{c ℓ} (K) : = A - B K C

. Let

C_{-}

denote the left-half plane, let

α > 0

and let

C_{α}

denote the set of all

z \in C

with

ℜ (z) \leq - α

, where

ℜ (z)

is the real part of z. Let

S^{q \times r}

denote the set of all matrices

K \in R^{q \times r}

, such that

A_{c ℓ}

is stable, i.e.,

σ (A_{c ℓ}) \subset C_{-}

(where

σ (A_{c ℓ})

is the spectrum of

A_{c ℓ}

). By

S_{α}^{q \times r}

, we denote the set of all matrices

K \in R^{q \times r}

, such that

σ (A_{c ℓ}) \subset C_{α}

. In this case, we say that

A_{c ℓ}

α

-stable. Below, we will occasionally write

S_{α}

instead of

S_{α}^{q \times r}

, when it is clear what the size of the related matrices is.

Let

K \in S_{α}^{q \times r}

be given. Substitution of

u = - K y = - K C x

into (2) gives

J (x_{0}, K) : = \int_{0}^{\infty} x {(t)}^{T} (Q + C^{T} K^{T} R K C) x (t) d t .

(3)

Since

Q + C^{T} K^{T} R K C > 0

and since

A_{c ℓ} (K)

is stable, it follows that the Lyapunov equation

A_{c ℓ} {(K)}^{T} P + P A_{c ℓ} (K) = - (Q + C^{T} K^{T} R K C)

(4)

has a unique solution

P > 0

, given by

P = - m a t ({(I_{p} \otimes A_{c ℓ} {(K)}^{T} + A_{c ℓ} {(K)}^{T} \otimes I_{p})}^{- 1} \cdot v e c (Q + C^{T} K^{T} R K C)),

(5)

where

v e c

puts all the columns of the given matrix in a single column and

m a t

is the inverse of

v e c

. Let us denote the solution (5) as

P (K)

. Substitution of (4) into (3) and noting that

\dot{x} (t) = A_{c ℓ} (K) x (t)

with

A_{c ℓ} (K)

stable, leads to

\begin{matrix} J (x_{0}, K) & = - \int_{0}^{\infty} x {(t)}^{T} (A_{c ℓ} {(K)}^{T} P (K) + P (K) A_{c ℓ} (K)) x (t) d t \\ = - \int_{0}^{\infty} \frac{d}{d t} (x {(t)}^{T} P (K) x (t)) d t = x_{0}^{T} P (K) x_{0} . \end{matrix}

Thus, we look for

K \in S_{α}^{q \times r}

that minimizes the functional

J (x_{0}, K) = x_{0}^{T} P (K) x_{0}

. When

x_{0}

is unknown, we seek

K \in S_{α}^{q \times r}

for which

σ_{m a x} (K) : = max (σ (P (K)))

(6)

is minimal. In this case, we get a robust LQR via SOF, in the sense that it minimizes

J (x_{0}, K)

for the worst possible (unknown)

x_{0}

. Note that

x_{0}^{T} P (K) x_{0} = {∥P {(K)}^{\frac{1}{2}} x_{0}∥}^{2} \leq {∥P {(K)}^{\frac{1}{2}}∥}^{2} {∥x_{0}∥}^{2} = ∥P (K)∥ {∥x_{0}∥}^{2} = σ_{m a x} (K) {∥x_{0}∥}^{2},

and that there exists

x_{0} \neq 0

for which equality holds. Therefore

\frac{J (x_{0}, K)}{{∥x_{0}∥}^{2}} \leq σ_{m a x} (K)

, where equality holds in the worst case.

Note that the functionals

J (x_{0}, K)

and

σ_{m a x} (K)

are generally not convex since their domain of definition

S^{q \times r}

(and therefore

S_{α}^{q \times r}

) is generally non-convex. Necessary conditions for optimality were given as three quadratic matrix equations in [2,3,4,5]. Necessary and sufficient conditions for optimality, based on linear matrix inequalities (LMI), were given in [6,7,8]. However, algorithms based on these formulations are generally not guaranteed to converge, seemingly because of the non-convexity of the coupled matrix equations or inequalities, and when they converge, it is to a local optimum only.

The application of SOFs in LQRs is appealing for several reasons: they are reliable and cheap, and their implementation is simple and direct. Moreover, the long-term memory of dynamic-feedbacks is useless for systems subject to random disturbances, to fast dynamic loadings or to impulses, and the application of state-feedbacks is not always possible, due to unavailability of full-state measurements (see [9], for example). On the other hand, in practical applications, the entries of the needed SOFs are bounded, and since the problem of SOFs with interval constrained entries is NP-hard (see [10,11]), one cannot expect the existence of a deterministic polynomial-time algorithm to solve this problem. Randomized algorithms are thus natural solutions to this problem. The probabilistic and randomized methods for the constrained SOF problem and robust stabilization via SOFs (among other hard problems) are discussed in [12,13,14,15]. The Ray-Shooting Method was recently introduced in [16], where it was utilized to derive the Ray-Shooting (RS) randomized algorithm for the minimal-gain SOF problem with regional pole-assignment, where the region can be non-convex and unconnected. For a survey of the SOF problem see [17] and for a recent survey of the robust SOF problem see [18].

The contribution of this research is as follows:

The suggested algorithm is based on the Ray-Shooting Method (see [16]), which, as opposed to smooth optimization methods, has the potential of finding a global optimum of continuous functions over compact non-convex and unconnected regions.
The suggested algorithm has a proof of convergence (in probability) and explicit complexity.
Experience with the algorithm shows good quality of controllers, high percent of success and good run-time for real-life systems. Thus, the suggested practical algorithm efficiently solves the problem of LQR via SOF.
The algorithm does not need to solve any Riccati equations and thus can be applied to large systems.
The suggested algorithm is one of the few that deals with LQR via SOF and has the ability to deal with discrete-time systems under the same formulation.

The reminder of the article is organized as follows:

In Section 2, we introduce the practical randomized algorithm for the problem of LQR via SOF. In Section 3, we give the results of the algorithm for some real-life systems and we compare its performance with the performance of a well known algorithm that has a proof of convergence to local minimum (under some reasonable assumptions). Finally, in Section 4 we conclude with some remarks.

2. The Practical Algorithm for the Problem of LQR via SOF

Assume that

K^{(0)} \in i n t (S_{α})

was found by the RS algorithm (see [16]) or by any other method (see [19,20,21]). Let

h > 0

and let

U^{(0)}

be a unit vector w.r.t. the Frobenius norm, i.e.,

{∥U^{(0)}∥}_{F} = 1

. Let

L^{(0)} = K^{(0)} + h \cdot U^{(0)}

and let

L

be the hyperplane defined by

L^{(0)} + V

, where

{〈V, U^{(0)}〉}_{F} = 0

. Let

r_{\infty} > 0

and let

R_{\infty}

denote the set of all

F \in L

, such that

{∥F - L^{(0)}∥}_{F} \leq r_{\infty}

. Let

R_{\infty} (ϵ) = R_{\infty} + \bar{B (0, ϵ)}

, where

\bar{B (0, ϵ)}

denotes the closed ball centered at 0 with radius

ϵ

(

0 < ϵ \leq \frac{1}{2}

), with respect to the Frobenius norm on

R^{q \times r}

. Let

D^{(0)} = C H (K^{(0)}, R_{\infty} (ϵ))

denote the convex-hull of the vertex

K^{(0)}

with the basis

R_{\infty} (ϵ)

. Let

S_{α}^{(0)} = S_{α} \cap D^{(0)}

and note that

S_{α}^{(0)}

is compact (but generally not convex). We wish to minimize the continuous function

σ_{m a x} (K)

(or the continuous function

J (x_{0}, K)

, when

x_{0}

is known) over the compact set

S_{α} \cap \bar{B (K^{(0)}, h)}

. Let

K_{*}

denote a point in

S_{α} \cap \bar{B (K^{(0)}, h)}

where the minimum of

σ_{m a x} (K)

is accepted. Obviously,

K_{*} \in D^{(0)}

, for some direction

U^{(0)}

, as above.

The suggested algorithm in Algorithm 1 works as follows:

We start with a point

K^{(0)} \in i n t (S_{α})

, found by the RS algorithm.

Assuming that

K_{*} \in D^{(0)}

, the inner-loop (

j = 1, \dots, n

) uses the Ray-Shooting Method in order to find an approximation of the global minimum of the function

σ_{m a x} (K)

over

S_{α}^{(0)}

—the portion of

S_{α}

bounded in the cone

D^{(0)}

. The proof of convergence in probability of the inner-loop and its complexity (under the above mentioned assumption) can be found in [16] (see also [22]). In the inner-loop, we choose a search direction by choosing a point F in

R_{\infty} (ϵ)

—the base of the cone

D^{(0)}

. Next, in the most inner-loop (

k = 1, \dots, s

) we scan the ray

K (t) : = (1 - t) K^{(0)} + t F

and record the best controller on it. Repeating this a sufficient number of times (as is given in (7) and in the discussion right after), we reach

K_{*}

(or an

ϵ

- neighborhood of it) with high probability, under the assumption that

K_{*} \in D^{(0)}

The outer-loop (

i = 1, \dots, m

) is used as a substitution for restarting the RS algorithm again and again, by taking

K^{(best)}

as the new vertex of the search cone instead of

K^{(0)}

and by choosing a different direction

U^{(0)}

. The choice of a different direction is made as a backup to the case where the above mentioned assumption didn’t hold in the previous iterations (see Remark 1 below). The replacement of

K^{(0)}

K^{(best)}

can be considered as a heuristic step, which is made instead of running the RS algorithm many times in order to generate “the best starting point”, which is relevant only if we actually evaluate

σ_{m a x} (K)

on each such point and take the point with the best value as the best starting point. Since we, in any case, evaluate

σ_{m a x} (K)

in the main algorithm, we could avoid the repeated execution of the RS algorithm. The outer-loop is similar to what is done in the Hide-And-Seek algorithm (see [23,24]). The convergence in probability of the Hide-And-Seek algorithm can be found in [25].

Remark 1.

The volume of

\bar{B (K^{(0)}, h)}

is given by

\frac{π^{ℓ / 2}}{Γ (ℓ / 2 + 1)} \cdot h^{ℓ}

where

ℓ : = q r

and Γ is the known Γ-function. The volume of

D^{(0)}

is given approximately (and exactly when

ϵ = 0

) by

\frac{h}{ℓ} \cdot \frac{π^{(ℓ - 1) / 2}}{Γ ((ℓ - 1) / 2 + 1)} \cdot r_{\infty}^{ℓ - 1}

. Thus, by taking

r_{\infty} = h

, the portion of

\bar{B (K^{(0)}, h)}

covered by

D^{(0)}

(i.e., the probability that

K_{*} \in D^{(0)}

) is given by

\frac{Γ (ℓ / 2 + 1)}{ℓ \cdot \sqrt{π} \cdot Γ ((ℓ - 1) / 2 + 1)}

. Let Θ denote the known relation between functions

f, g : N \to R

defined by

f (n) = Θ (g (n))

if and only if

{lim}_{n \to \infty} \frac{f (n)}{g (n)} = 1

. Since

Γ (ℓ / 2 + 1) = Θ (\sqrt{2 π} e^{- ℓ / 2} {(\frac{ℓ}{2})}^{\frac{ℓ + 1}{2}})

and since

{(\frac{ℓ - 1}{ℓ})}^{ℓ / 2} \to e^{- 1 / 2}

when

ℓ \to + \infty

, it follows that

\frac{Γ (ℓ / 2 + 1)}{ℓ \cdot \sqrt{π} \cdot Γ ((ℓ - 1) / 2 + 1)} = Θ (\frac{1}{e \cdot \sqrt{2 π ℓ}}) .

Therefore, by taking

m = ⌈e \cdot \sqrt{2 π ℓ}⌉

iterations in the outer-loop, we have

K_{*} \in D^{(0)}

almost surely. Specifically, when

ℓ \geq 12

, we suggest taking

m = 2 ℓ

and

ℓ \times ℓ

orthogonal matrix

U = [\begin{matrix} u_{1} & u_{2} & \dots & u_{ℓ} \end{matrix}]

, and to take the directions

U_{j}^{(0)} = \pm m a t (u_{j}), j = 1, \dots, ℓ

in the outer-loop.

The complexity of the suggested practical algorithm measured as the number of its arithmetic operations is given as follows:

computing the matrix $P (K (t))$ as in (5) takes $O ({(p^{2})}^{3}) = O (p^{6})$ , since the dominant operation is the inversion of the $p^{2} \times p^{2}$ matrix there.
checking $K (t) \in S_{α}$ by checking the $α$ -stability of $A_{c l} (K (t))$ (as well as computing $σ (K (t))$ ), takes $O (p^{3})$ for computing the characteristic polynomial of $A_{c ℓ} (K (t))$ (of $P (K (t))$ , respectively) and $O (p {log}_{2}^{2} (p) ({log}_{2}^{2} (p) + {log}_{2}^{2} (b)))$ for computing approximations ${\tilde{λ}}_{j}$ for all the eigenvalues $λ_{j}, j = 1, \dots, p$ of $A_{c ℓ} (K (t))$ (of $P (K (t))$ , respectively), with accuracy $|{\tilde{λ}}_{j} - λ_{j}| < 2^{2 - \frac{b}{p}}$ where $b \geq p {log}_{2} (p)$ . The approximated eigenvalues can be computed to the accuracy $2^{2 - \frac{b}{p}} = ϵ$ , with $b = (2 + {log}_{2} (\frac{1}{ϵ})) p$ , by the algorithm of V. Y. Pan (see [26]). We end up with $O (p^{3})$ for these operations.
computing uniformly distributed $q \times r$ matrix takes $O (max {(q, r)}^{3})$ operations.

We therefore have a total complexity of

O (m n s (max {(q, r)}^{3} + p^{6}))

Let a closed

ϵ

-neighborhood of

K_{*}

D^{(0)}

be defined by

S_{α}^{(0)} (ϵ) = \{K \in S_{α}^{(0)} |σ (K) \leq σ (K_{*}) + ϵ\} .

Let the idealized algorithm be the algorithm that samples the search space

D^{(0)}

until hitting

S_{α}^{(0)} (ϵ)

, where the sampling is according to a general p.d.f. g and a related generator G. For

0 < β < 1

, the number of iterations needed to guarantee a probability of at least

1 - β

to hit

S_{α}^{(0)} (ϵ)

is given by

⌈\frac{M_{g} V o ℓ (D^{(0)})}{m_{g} V o ℓ (S_{α}^{(0)} (ϵ))} |ln (β)|⌉,

(7)

where

Vo ℓ

denotes the volume of the related set and

M_{g}, m_{g}

are the essential-supremum and essential-infimum of g over

D^{(0)}

, respectively (see [16]). Similarly to what is done in [16], one can show that the last is

O (\frac{|ln (β)| h M_{g}}{ϵ m_{g}} {(\frac{r_{\infty}}{r_{ϵ}})}^{q r})

, where

r_{ϵ}

is a radius of a ball of a basis of a cone with height

ϵ

and vertex

K_{*}

that has a volume that equals to

V o ℓ (S_{α}^{(0)} (ϵ))

. This results in an exponential number of iterations, but if we restrict the input of the algorithm to systems with

q, r

satisfying

q \leq q_{0}, r \leq r_{0}

when

q_{0}, r_{0}

are fixed, then, the number of iterations would be

O (\frac{|ln (β)| h M_{g}}{ϵ m_{g}} {(\frac{r_{\infty}}{r_{ϵ}})}^{q_{0} r_{0}})

, i.e., polynomial in

(\frac{r_{\infty}}{r_{ϵ}})

—which can be considered as the true size of the problem (for a fixed

p, q, r

). In this sense we can say that the algorithm is efficient. The total number of arithmetic operations of the idealized algorithm that guarantees a probability at least

1 - β

to hit

S_{α}^{(0)} (ϵ)

is therefore given by

O (\frac{|ln (β)| h}{ϵ} {(\frac{r_{\infty}}{r_{ϵ}})}^{q_{0} r_{0}} (max {(q, r)}^{3} + p^{6}))

, since sampling points according to the uniform distribution g (and therefore

m_{g} = M_{g} = 1

) and the related generator G, takes

O (max {(q, r)}^{3})

For the sake of comparison, which will be presented in the next section, we bring here, in Algorithm 2, the algorithm of D. Moerder and A. Calise (see [5]) adjusted to our formulation of the problem, which we call: the MC Algorithm. To the best of our knowledge, this is the best algorithm for LQR via SOF published so far.

Algorithm 1: The practical randomized algorithm for the LQR via static-output-feedback (SOF) problem.

Input:

0 < ϵ \leq \frac{1}{2}, α, h, r_{\infty} > 0

, integers:

m, n, s > 0

controllable pairs

(A, B)

and

(A^{T}, C^{T})

matrices

Q > 0, R \geq 0

and

K^{(0)} \in i n t (S_{α})

Output:

K \in S_{α}

close as possible to

K_{*}

1. compute

P (K^{(0)})

as in (5)

P^{(best)} \leftarrow P (K^{(0)})

σ_{m a x}^{(best)} \leftarrow max (σ (P^{(best)}))

4. for

i = 1

to m do

4.1. choose

U^{(0)}

such that

{∥U^{(0)}∥}_{F} = 1

uniformly at random

4.2. let

L^{(0)} \leftarrow K^{(0)} + h \cdot U^{(0)}

4.3. for

j = 1

to n do

4.3.1. choose

F \in R_{\infty} (ϵ)

uniformly at random

4.3.1.1. for

k = 1

to s do

4.3.1.1.1.

t \leftarrow \frac{k}{s}

4.3.1.1.2.

K (t) \leftarrow (1 - t) K^{(0)} + t F

4.3.1.1.3. if

K (t) \in S_{α}

then

4.3.1.1.3.1. compute

P (K (t))

as in (5)

4.3.1.1.3.2.

σ_{m a x} (K (t)) \leftarrow max (σ (P (K (t))))

4.3.1.1.3.3. if

(σ_{m a x} (K (t)) < σ_{m a x}^{(best)})

then

4.3.1.1.3.3.1.

K^{(best)} \leftarrow K (t)

4.3.1.1.3.3.2.

P^{(best)} \leftarrow P (K (t))

4.3.1.1.3.3.3.

σ_{m a x}^{(best)} \leftarrow σ_{m a x} (K (t))

4.4.

K^{(0)} \leftarrow K^{(best)}

5. return

K^{(best)}, P^{(best)}, σ_{m a x}^{(best)}

Algorithm 2: The MC Algorithm.

Input:

0 < ϵ, α

, integers:

m, s > 0

controllable pairs

(A, B)

and

(A^{T}, C^{T})

matrices

Q > 0, R > 0

and

K_{0} \in i n t (S_{α})

Output:

K \in S_{α}

close as possible to

K_{*}

j \leftarrow 0

A_{0} \leftarrow A - B K_{0} C

P_{0} \leftarrow - m a t (\begin{matrix} {(I_{p} \otimes A_{0}^{T} + A_{0}^{T} \otimes I_{p})}^{- 1} \cdot \\ v e c (Q + C^{T} K_{0}^{T} R K_{0} C) \end{matrix})

S_{0} \leftarrow - m a t ({(I_{p} \otimes A_{0} + A_{0} \otimes I_{p})}^{- 1} \cdot v e c (I_{p}))

σ_{m a x} (K_{0}) \leftarrow max (σ (P_{0}))

Δ K_{0} \leftarrow R^{- 1} B^{T} P_{0} S_{0} C^{T} {(C S_{0} C)}^{- 1} - K_{0}

f l a g \leftarrow 0

8. for

k = 1

to s do

8.1.

t \leftarrow \frac{k}{s}

8.2.

K (t) \leftarrow (1 - t) K_{0} + t Δ K_{0}

8.3. if

K (t) \in S_{α}

then

8.3.1.

A (t) \leftarrow A - B K (t) C

8.3.2.

P (t) \leftarrow - m a t (\begin{matrix} {(I_{p} \otimes A {(t)}^{T} + A {(t)}^{T} \otimes I_{p})}^{- 1} \cdot \\ v e c (Q + C^{T} K {(t)}^{T} R K (t) C) \end{matrix})

8.3.3.

S (t) \leftarrow - m a t ({(I_{p} \otimes A (t) + A (t) \otimes I_{p})}^{- 1} \cdot v e c (I_{p}))

8.3.4.

σ_{m a x} (K (t)) \leftarrow max (σ (P (t)))

8.3.5. if

σ_{m a x} (K (t)) < σ_{m a x} (K_{0})

then

8.3.5.1.

K_{1} \leftarrow K (t)

8.3.5.2.

A_{1} \leftarrow A - B K_{1} C

8.3.5.3.

P_{1} \leftarrow P (t)

8.3.5.4.

S_{1} \leftarrow S (t)

8.3.5.5.

σ_{m a x} (K_{1}) \leftarrow σ_{m a x} (K (t))

8.3.5.6.

f l a g \leftarrow 1

9. if

f l a g = = 1

then

9.1. while

(|σ_{m a x} (K_{j + 1}) - σ_{m a x} (K_{j})| \geq ϵ)

and

(j < m)

9.1.1.

Δ K_{j} \leftarrow R^{- 1} B^{T} P_{j} S_{j} C^{T} {(C S_{j} C)}^{- 1} - K_{j}

9.1.2. for

k = 1

to s do

9.1.2.1.

t \leftarrow \frac{k}{s}

9.1.2.2.

K (t) \leftarrow (1 - t) K_{j} + t Δ K_{j}

9.1.2.3. if

K (t) \in S_{α}

then

9.1.2.3.1.

A (t) \leftarrow A - B K (t) C

9.1.2.3.2.

P (t) \leftarrow - m a t (\begin{matrix} {(I_{p} \otimes A {(t)}^{T} + A {(t)}^{T} \otimes I_{p})}^{- 1} \cdot \\ v e c (Q + C^{T} K {(t)}^{T} R K (t) C) \end{matrix})

9.1.2.3.3.

S (t) \leftarrow - m a t ({(I_{p} \otimes A (t) + A (t) \otimes I_{p})}^{- 1} \cdot v e c (I_{p}))

9.1.2.3.4.

σ_{m a x} (K (t)) \leftarrow max (σ (P (t)))

9.1.2.3.5. if

σ_{m a x} (K (t)) < σ_{m a x} (K_{j})

then

9.1.2.3.5.1.

K_{j + 1} \leftarrow K (t)

9.1.2.3.5.2.

A_{j + 1} \leftarrow A - B K_{j + 1} C

9.1.2.3.5.3.

P_{j + 1} \leftarrow P (t)

9.1.2.3.5.4.

S_{j + 1} \leftarrow S (t)

9.1.2.3.5.5.

σ_{m a x} (K_{j + 1}) \leftarrow σ_{m a x} (K (t))

9.1.2.3.5.6.

j \leftarrow j + 1

10. return

K^{(best)} \leftarrow K_{j}, A_{j}, P_{j}, S_{j}, σ_{m a x}^{(best)} \leftarrow σ_{m a x} (K^{(b e s t)})

3. Experiments

In the following experiments we applied the Algorithm 1 and Algorithm 2, on systems taken from the liberaries [27,28,29]. We took only the systems with controllable

(A, B), (A^{T}, C^{T})

pairs, for which the RS algorithm succeeded in finding SOFs (see [16], Table 8, p. 231). In order to initialize the MC Algorithm, we also used the RS algorithm to find a starting

α

-stabilizing static-output-feedback, for known optimal value for

α

. In all the experiments for the suggested algorithm we used

m = 100

n = 100

s = 100

h = 100, r_{\infty} = 100, ϵ = 10^{- 16}

, and for the MC Algorithm we used

m = 10, 000, s = 100

(in order to get the same number of

10^{6}

overall iterations and the same number

s = 100

of iterations for the local search). In every case, we took

Q = I_{p}, R = I_{q}

. The Stability Margin column of Table 1 relates to

α > 0

for which the real part of any eigenvalue of the closed-loop is less than or equal to

- α

. The values of

α

in Table 1 relates to the largest

α

for which the RS algorithm succeeded in finding

K^{(0)}

. The reason is that if

λ

is an eigenvalue of

A_{c ℓ} (K)

with corresponding eigenvector v then, (4) implies

ℜ (λ) = - \frac{v^{*} (Q + C^{T} K^{T} R K C) v}{2 v^{*} P (K) v} \leq - \frac{v^{*} Q v}{2 v^{*} P (K) v} .

It follows that minimizing

σ_{m a x} (K)

results in a larger abscissa. Thus, it is worth searching for a starting point

K^{(0)}

that maximizes the abscissa

α

. This can be done efficiently by running a binary search on the abscissa and using the RS algorithm as an oracle. Note that RS CPU time appearing in the third column of Table 1 relates to running the RS algorithm for known optimal value of the abscissa. The RS algorithm is sufficiently fast also for this purpose, but other methods such as the HIFOO algorithm (see [19]) can be applied for this purpose as well.

Let

σ_{m a x} (F)

denote the functional (6) for the system

(A, B, I_{p})

, where

A - B F

is stable, i.e.,

F \in S^{q \times p}

. Let

P (F)

denote the unique solution of (4) for the system

(A, B, I_{p})

with F as above. Let

σ_{m a x} (K)

denote the functional (4) for the system

(A, B, C)

with

K \in S^{q \times r}

and related Lyapunov matrix

P = P (K)

. Now, if

A - B K C

is stable for some K then,

A - B F

is stable for

F = K C

(but there might exist a F such that

A - B F

is stable, but which cannot be defined as

K C

for some

q \times r

matrix K). Therefore

σ_{m a x} (F_{*}) = min_{F \in S^{q \times p}} σ_{m a x} (F) \leq min_{K \in S_{α}^{q \times r} \cap \bar{B (K^{(0)}, h)}} σ_{m a x} (K) = σ_{m a x} (K_{*}),

(8)

where

F_{*}

is an optimal LQR state-feedback controller for the system

(A, B, I_{p})

. We conclude that

σ_{m a x} (F_{*}) \leq σ_{m a x} (K_{*}) \leq σ_{m a x} (K^{(best)})

. Thus,

σ_{m a x} (F_{*})

is a lower-bound for

σ_{m a x} (K^{(best)})

and can serve as a good estimator for it (as is evidently seen from Table 1 in many cases) in order to stop the algorithm earlier, since

σ_{m a x} (F_{*})

can be calculated in advance.

The fourth column of Table 1 represents

σ_{m a x} (K^{(0)})

. The fifth column stands for

σ_{m a x} (K^{(best)})

, where the number in the parentheses is the relative improvement over

σ_{m a x} (K^{(0)})

, in percent. The sixth column is for

σ_{m a x} (F_{*})

and the seventh column is for the CPU time of the suggested algorithm, given in seconds. The eighth column stands for

σ_{m a x} (K^{(best)})

found by the MC Algorithm and finally, the ninth column stands for the CPU time of the MC Algorithm.

Regarding the suggested algorithm, note that the relative improvement of

σ_{m a x}^{(best)}

over

σ_{m a x}^{(0)}

is at least

1000 %

for the systems: AC6, AC11, AC12, DIS4, REA1, TF1, NN9, NN15 and NN17 (i.e., in 9 out of 29 systems). The AC12 should be noted, for which the improvement is

7.55 \cdot 10^{17} %

! Note also that for the (13 out of 29) systems: AC1, AC2, AC6, AC11, AC12, AC15, HE3, DIS4, REA1, REA2, HF2D10, HF2D11 and NN15, the value of

σ_{m a x}^{(best)}

for

(A, B, C)

is very close to the value of

σ_{m a x} (F_{*})

for

(A, B)

Regarding the comparison with the MC Algorithm, we conclude that the MC algorithm actually failed in finding any improvement of

σ_{m a x}^{(best)}

over

σ_{m a x}^{(0)}

for the systems: AC11, HE4, ROC1, ROC2, TF1 and NN12. The suggested algorithm performs significantly better regarding this key of performance for the systems: AC6, AC11, AC15, HE3, HE4, TF1, ROC4, NN9 and NN12. The suggested algorithm performs slightly better for the systems: AC1, AC2, DIS4, DIS5, HE1, HF2D10, HF2D11, REA1, REA2, ROC1, TMD, NN1, NN5 and NN13. Finally, the MC algorithm performs slightly better only for the systems: AC5, AC12 and NN16. We conclude that the suggested algorithm outperforms the MC Algorithm regarding the above-mentioned key of performance, although the MC Algorithm outperforms the suggested algorithm in terms of CPU time.

The MC Algorithm seems to perform better locally, while the suggested algorithm seems to perform better globally. Thus, practically, the best approach would be to apply the suggested algorithm in order to find a close neighborhood of a global minimum and then to apply the MC Algorithm on the result, for the local optimization.

4. Concluding Remarks

For a discrete-time system

\{\begin{matrix} x_{k + 1} = A x_{k} + B u_{k} \\ y_{k} = C x_{k} \end{matrix}

(9)

and cost functional

J (x_{0}, u) : = \sum_{k = 0}^{\infty} (x_{k}^{T} Q x_{k} + u_{k}^{T} R x_{k}),

(10)

let

u_{k} = - K y_{k}

be the SOF, and let

A_{c ℓ} (K) : = A - B K C

be the closed-loop matrix. Let

D

denote the open unit-disk, let

0 < α < 1

and let

D_{α}

denote the set of all

z \in D

with

|z| \leq 1 - α

(where

|z|

is the absolute value of z). Let

S^{q \times r}

denote the set of all matrices

K \in R^{q \times r}

such that

σ (A_{c ℓ}) \subset D

(i.e., stable in the discrete-time sense), and let

S_{α}^{q \times r}

denote the set of all matrices

K \in R^{q \times r}

such that

σ (A_{c ℓ}) \subset D_{α}

. In this case we say that

A_{c ℓ}

α

-stable. Let

K \in S_{α}^{q \times r}

be given. Substitution of

u_{k} = - K y_{k} = - K C x_{k}

into (10) yields

J (x_{0}, K) : = \sum_{k = 0}^{\infty} x_{k}^{T} (Q + C^{T} K^{T} R K C) x_{k} .

(11)

Since

Q + C^{T} K^{T} R K C > 0

and since

A_{c ℓ} (K)

is stable, it follows that the Stein equation

P - A_{c ℓ} {(K)}^{T} P A_{c ℓ} (K) = Q + C^{T} K^{T} R K C

(12)

has a unique solution

P > 0

, given by

P (K) = m a t ({(I_{p} \otimes I_{p} - A_{c ℓ} {(K)}^{T} \otimes A_{c ℓ} {(K)}^{T})}^{- 1} \cdot v e c (Q + C^{T} K^{T} R K C)) .

Substitution of (12) into (11) and noting that

x_{k} = A_{c ℓ} {(K)}^{k} x_{0}

with

A_{c ℓ} (K)

stable, leads to

\begin{matrix} J (x_{0}, K) & = \sum_{k = 0}^{\infty} x_{k}^{T} (P - A_{c ℓ} {(K)}^{T} P A_{c ℓ} (K)) x_{k} \\ = \sum_{k = 0}^{\infty} x_{0}^{T} A_{c ℓ} {(K)}^{T k} (P - A_{c ℓ} {(K)}^{T} P A_{c ℓ} (K)) A_{c ℓ} {(K)}^{k} x_{0} = x_{0}^{T} P (K) x_{0} . \end{matrix}

Thus, we look for

K \in S_{α}^{q \times r}

that minimizes the functional

J (x_{0}, K) = x_{0}^{T} P (K) x_{0}

, and when

x_{0}

is unknown, we seek

K \in S_{α}^{q \times r}

for which

σ_{m a x} (K) : = max (σ (P (K)))

is minimal. Similarly to the continuous-time case, we have

\frac{J (x_{0}, K)}{{∥x_{0}∥}^{2}} \leq σ_{m a x} (K)

with equality in the worst case. Finally, if

λ

is an eigenvalue of

A_{c ℓ} (K)

and v is a corresponding eigenvector then (12) yields

1 - {|λ|}^{2} = \frac{v^{*} (Q + C^{T} K^{T} R K C) v}{v^{*} P (K) v} \geq \frac{v^{*} Q v}{v^{*} P (K) v}

. Therefore

{|λ|}^{2} \leq 1 - \frac{v^{*} Q v}{v^{*} P (K) v}

, and thus, minimizing

σ_{m a x} (K)

results in larger abscissa. Now, the suggested algorithm can be readily applied to discrete-time systems. As to the MC Algorithm, we are not aware of any discrete-time analogue of it.

We conclude that the Ray-Shooting Method is a powerful tool, since it practically solves the problem of LQR via SOF, for real-life systems. The suggested algorithm has good performance, and is proved to converge in probability. The suggested method can similarly handle the problem of LQR via SOF for discrete-time systems. Obviously, this enlarges the scope and usability of the Ray-Shooting Method and we expect to receive more results in this direction.

Acknowledgments

I wish to dedicate this article to the memory of my beloved father Pinhas Peretz and to the memory of my beloved wife and friend Rivka Rimonda Peretz.

Conflicts of Interest

The author declares no conflict of interest.

References

Peretz, Y. A characterization of all the static stabilizing controllers for LTI systems. Linear Algebra Its Appl. 2012, 437, 525–548. [Google Scholar] [CrossRef]
Johnson, T.; Atahns, M. On the design of optimal dynamic compansators for linear constant systems. IEEE Trans. Autom. Control 1970, 15, 658–660. [Google Scholar] [CrossRef]
Levine, W.; Athans, M. On the determination of the optimal constant output feedback gains for linear multivariables systems. IEEE Trans. Autom. Control 1970, 15, 44–48. [Google Scholar] [CrossRef]
Levine, W.; Johnson, T.L.; Athans, M. Optimal limited state variable feedback controllers for linear systems. IEEE Trans. Autom. Control 1971, 16, 785–793. [Google Scholar] [CrossRef]
Moerder, D.; Calise, A. Convergence of numerical algorithm for calculating optimal output feedback gains. IEEE Trans. Autom. Control 1985, 30, 900–903. [Google Scholar] [CrossRef]
Iwasaki, T.; Skelton, R. All controllers for the general H_∞ control problem: LMI existance conditions and stste space formulas. Automatica 1994, 30, 1307–1317. [Google Scholar] [CrossRef]
Iwasaki, T.; Skelton, R.E. Linear quadratic suboptimal control with static output feedback. Syst. Control Lett. 1994, 23, 421–430. [Google Scholar] [CrossRef]
Peres, P.L.D.; Geromel, J.; de Souza, S. Optimal H₂ control by output feedback. In Proceedings of the 32nd IEEE Conference on Decision and Control, San Antonio, TX, USA, 15–17 December 1993; pp. 102–107. [Google Scholar]
Camino, J.F.; Zampieri, D.E.; Peres, P.L.D. Design of A Vehicular Suspension Controller by Static Output Feedback. In Proceedings of the American Control Conference, San Diego, CA, USA, 2–4 June 1999. [Google Scholar]
Blondel, V.; Tsitsiklis, J.N. NP-hardness of some linear control design problems. SIAM J. Control Optim. 1997, 35, 2118–2127. [Google Scholar] [CrossRef]
Nemirovskii, A. Several NP-hard problems arising in robust stability analysis. Math. Control Signals Syst. 1993, 6, 99–105. [Google Scholar] [CrossRef]
Arzelier, D.; Gryazina, E.N.; Peaucelle, D.; Polyak, B.T. Mixed LMI/randomized methods for static output feedback control. In Proceedings of the American Control Conference, Baltimore, MD, USA, 30 June–2 July 2010; pp. 4683–4688. [Google Scholar]
Tempo, R.; Calafiore, G.; Dabbene, F. Randomized Algorithms for Analysis and Control of Uncertain Systems; Springer: London, UK, 2005. [Google Scholar]
Tempo, R.; Ishii, H. Monte Carlo and Las Vegas Randomized Algorithms for Systems and Control. Eur. J. Control 2007, 13, 189–203. [Google Scholar] [CrossRef]
Vidyasagar, M.; Blondel, V.D. Probabilistic solutions to some NP-hard matrix problems. Automatica 2001, 37, 1397–1405. [Google Scholar] [CrossRef]
Peretz, Y. A randomized approximation algorithm for the minimal-norm static-output-feedback problem. Automatica 2016, 63, 221–234. [Google Scholar] [CrossRef]
Syrmos, V.L.; Abdallah, C.; Dorato, P.; Grigoradis, K. Static Output Feedback: A Survey. Automatica 1997, 33, 125–137. [Google Scholar] [CrossRef]
Sadabadi, M.S.; Peaucelle, D. From static output feedback to structured robust static output feedback: A survey. Annu. Rev. Control 2016, 42, 11–26. [Google Scholar] [CrossRef]
Gumussoy, S.; Henrion, D.; Millstone, M.; Overton, M.L. Multiobjective Robust Control with HIFOO 2.0. In Proceedings of the IFAC Symposium on Robust Control Design, Haifa, Israel, 16–18 June 2009. [Google Scholar]
Yang, K.; Orsi, R. Generalized pole placement via static output feedback: A methodology based on projections. Automatica 2006, 42, 2143–2150. [Google Scholar] [CrossRef]
Henrion, D.; Loefberg, J.; Kočvara, M.; Stingl, M. Solving Polynomial static output feedback problems with PENBMI. In Proceedings of the IEEE Conference on Decision and Control, Sevilla, Spain, 15 December 2005. [Google Scholar]
Peretz, Y. On applications of the Ray-Shooting method for structured and structured-sparse static-output-feedbacks. Int. J. Syst. Sci. 2017, 48, 1902–1913. [Google Scholar]
Zabinsky, Z.B. Stochastic Adaptive Search for Global Optimization; Kluer Academic Publishers: Dordrecht, The Netherlands, 2003. [Google Scholar]
Romeijn, H.E.; Smith, R.L. Simulated Annealing for Constrained Global Optimization. J. Glob. Optim. 1994, 5, 101–126. [Google Scholar] [CrossRef]
Bélisle, C.J.P. Convergence Theorems for a Class of Simulated Annealing Algorithms on $R^{d}$ . J. Appl. Probab. 1992, 29, 885–895. [Google Scholar]
Pan, V.Y. Univariate polynomials: Nearly optimal algorithms for numerical factorization and root-finding. J. Symb. Comput. 2002, 33, 701–733. [Google Scholar] [CrossRef]
Leibfritz, F. COMPl_eib: Constrained Matrix-Optimization Problem Library—A Collection of Test Examples for Nonlinear Semidefinite Programs, Control System Design and Related Problems; Technical Report; Department of Methematics, University of Trier: Trier, Germany, 2003. [Google Scholar]
Leibfritz, F. Description of the Benchmark Examples in COMPl_eib 1.0; Technical Report; Department of Methematics, University of Trier: Trier, Germany, 2003. [Google Scholar]
Leibfritz, F.; Lipinski, W. COMPl_eib 1.0—User Manual and Quick Reference; Technical Report; Department of Methematics, University of Trier: Trier, Germany, 2004. [Google Scholar]

Table 1. Results of the suggested randomized algorithm for LQR via SOF compared with the results of the MC Algorithm.

System	Stability Margin	RS CPU Time	$σ_{\max}^{(0)}$ for $(A, B, C)$	$σ_{\max}^{(best)}$ for $(A, B, C)$ Suggested Algorithm	$σ_{\max} (F_{*})$ for $(A, B)$	Suggested Algorithm CPU Time	$σ_{\max}^{(best)}$ for $(A, B, C)$ MC Algorithm	MC Algorithm CPU Time
AC1	$0.1$	0.0625	20.9727	$\begin{matrix} 13.3987 \\ (56.52 %) \end{matrix}$	13.0686	43.1718	16.2315	0.1562
AC2	$0.1$	0.0312	20.9727	$\begin{matrix} 13.4217 \\ (56.25 %) \end{matrix}$	13.0686	42.3906	16.2315	0.2500
AC5	$0.875$	0.0625	$2.2821 \times 10^{6}$	$\begin{matrix} 2.2208 \times 10^{6} \\ (2.76 %) \end{matrix}$	$8.4264 \times 10^{5}$	56.8125	$2.0608 \times 10^{6}$	0.1093
AC6	$0.875$	0.1093	463.8502	$\begin{matrix} 6.2413 \\ (7.33 \times 10^{3} %) \end{matrix}$	5.9721	43.8125	91.3276	0.0625
AC11	$0.1$	0.1093	$1.0661 \times 10^{3}$	$\begin{matrix} 7.4244 \\ (1.42 \times 10^{4} %) \end{matrix}$	5.8648	31.6406	$1.0661 \times 10^{3}$	0.0625
AC12	$0.1$	< $10^{- 4}$	$4.1640 \times 10^{19}$	$\begin{matrix} 5.5149 \times 10^{3} \\ (7.55 \times 10^{17} %) \end{matrix}$	$2.7690 \times 10^{3}$	46.3437	$2.9950 \times 10^{3}$	3.0625
AC15	$0.25$	< $10^{- 4}$	112.0500	$\begin{matrix} 106.1312 \\ (5.57 %) \end{matrix}$	104.8458	42.0000	111.9781	0.1562
HE1	$0.1$	< $10^{- 4}$	26.7550	$\begin{matrix} 9.1169 \\ (1.93 \times 10^{2} %) \end{matrix}$	2.9961	44.8125	12.9412	0.1250
HE3	$0.25$	0.0312	944.3921	$\begin{matrix} 632.1294 \\ (49.39 %) \end{matrix}$	611.8468	100.7812	668.0894	0.4218
HE4	$0.01$	0.0625	$4.7440 \times 10^{3}$	$\begin{matrix} 474.7121 \\ (8.99 \times 10^{2} %) \end{matrix}$	229.9171	152.9531	$4.7440 \times 10^{3}$	< $10^{- 4}$
ROC1	$10^{- 16}$	0.0468	$1.4211 \times 10^{4}$	$\begin{matrix} 1.1368 \times 10^{4} \\ (25 %) \end{matrix}$	$1.1207 \times 10^{3}$	81.2500	$1.4211 \times 10^{4}$	0.0937
ROC4	$10^{- 16}$	0.1875	$1.7688 \times 10^{4}$	$\begin{matrix} 8.6513 \times 10^{3} \\ (1.04 \times 10^{2} %) \end{matrix}$	$8.5454 \times 10^{2}$	79.7968	$1.7688 \times 10^{4}$	0.0468
DIS4	1	< $10^{- 4}$	$3.2208 \times 10^{7}$	$\begin{matrix} 1.7529 \\ (1.83 \times 10^{9} %) \end{matrix}$	1.7504	68.2812	2.0166	0.7031
DIS5	$0.1$	0.0312	$3.9985 \times 10^{5}$	$\begin{matrix} 2.3800 \times 10^{5} \\ (68 %) \end{matrix}$	$9.0756 \times 10^{4}$	$109.8906$	$2.8304 \times 10^{5}$	0.3593
REA1	$0.75$	0.0468	$9.3478 \times 10^{5}$	$\begin{matrix} 2.2790 \\ (4.10 \times 10^{7} %) \end{matrix}$	2.2265	58.3906	2.7385	0.5156
REA2	$0.1$	0.0312	9.8349	$\begin{matrix} 2.2640 \\ (3.34 \times 10^{2} %) \end{matrix}$	2.2443	69.9843	2.7770	0.3125
TMD	$0.05$	0.0312	45.0958	$\begin{matrix} 27.2080 \\ (65.74 %) \end{matrix}$	16.7680	37.2031	27.6023	0.8281
TF1	$10^{- 16}$	0.1093	$9.1632 \times 10^{3}$	$\begin{matrix} 193.2880 \\ (4.64 \times 10^{3} %) \end{matrix}$	58.1296	51.3750	$9.1632 \times 10^{3}$	0.0312
HF2D10	$10^{- 16}$	0.0312	1.8090	$\begin{matrix} 1.4032 \\ (28.91 %) \end{matrix}$	1.3832	118.2656	1.4269	2.0781
HF2D11	$10^{- 16}$	0.0312	0.4315	$\begin{matrix} 0.3699 \\ (16.65 %) \end{matrix}$	0.3676	154.0312	0.3784	1
NN1	$10^{- 16}$	0.0312	$7.4631 \times 10^{3}$	$\begin{matrix} 1.6144 \times 10^{3} \\ (3.62 \times 10^{2} %) \end{matrix}$	106.7801	50.4375	$2.3561 \times 10^{3}$	0.2500
NN5	$0.01$	0.0468	$3.9123 \times 10^{4}$	$\begin{matrix} 9.6741 \times 10^{3} \\ (3.04 \times 10^{2} %) \end{matrix}$	$2.8787 \times 10^{3}$	129.2812	$9.8102 \times 10^{3}$	0.3125
NN9	$0.01$	0.0312	$4.0577 \times 10^{3}$	$\begin{matrix} 295.6797 \\ (1.27 \times 10^{3} %) \end{matrix}$	21.2937	40.6093	$3.7349 \times 10^{3}$	0.0937
NN12	$0.01$	0.0937	934.5127	$\begin{matrix} 236.4467 \\ (2.95 \times 10^{2} %) \end{matrix}$	30.3714	35.3281	934.5127	0.0625
NN13	$0.1$	0.0312	6.4782	$\begin{matrix} 1.7094 \\ (2.66 \times 10^{2} %) \end{matrix}$	0.6299	33.9218	1.8055	0.4531
NN14	$0.1$	0.0312	4.8907	$\begin{matrix} 1.6664 \\ (1.93 \times 10^{2} %) \end{matrix}$	0.6299	33.4687	1.7922	0.1718
NN15	$0.01$	< $10^{- 4}$	$1.3309 \times 10^{5}$	$\begin{matrix} 387.3778 \\ (3.42 \times 10^{4} %) \end{matrix}$	386.5741	94.2031	387.4183	0.3281
NN16	$0.1$	0.1093	35.9487	$\begin{matrix} 6.0864 \\ (4.90 \times 10^{2} %) \end{matrix}$	2.3276	45.0468	5.9982	0.7343
NN17	$0.1$	0.0312	$2.6404 \times 10^{3}$	$\begin{matrix} 36.7664 \\ (7.08 \times 10^{3} %) \end{matrix}$	3.1308	28.3281	666.7218	0.0937

© 2018 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Peretz, Y. On Application of the Ray-Shooting Method for LQR via Static-Output-Feedback. Algorithms 2018, 11, 8. https://doi.org/10.3390/a11010008

AMA Style

Peretz Y. On Application of the Ray-Shooting Method for LQR via Static-Output-Feedback. Algorithms. 2018; 11(1):8. https://doi.org/10.3390/a11010008

Chicago/Turabian Style

Peretz, Yossi. 2018. "On Application of the Ray-Shooting Method for LQR via Static-Output-Feedback" Algorithms 11, no. 1: 8. https://doi.org/10.3390/a11010008

APA Style

Peretz, Y. (2018). On Application of the Ray-Shooting Method for LQR via Static-Output-Feedback. Algorithms, 11(1), 8. https://doi.org/10.3390/a11010008

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

On Application of the Ray-Shooting Method for LQR via Static-Output-Feedback

Abstract

1. Introduction

2. The Practical Algorithm for the Problem of LQR via SOF

3. Experiments

4. Concluding Remarks

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI