**2. Preliminaries**

The Ray-Shooting Method was recently introduced in [10], where it was used

• LQR optimal, *H*∞-optimal, and *H*2-optimal proportional-integral-differential

The contribution of the research presented in the current chapter is as follows:

1.The randomized algorithm presented here (which we call the RS algorithm) is based on the Ray-Shooting Method (see [10]), which opposed to smooth optimization methods, and has the potential of finding a global optimum of continuous functions over compact non-convex and unconnected regions.

2.The RS algorithm has a proof of convergence in probability and explicit

4.The RS algorithm does not need to solve any Riccati or quadratic matrix

5.The RS algorithm is one of the few, dealing with the problem of LQR via SOF

6.A deterministic algorithm for the problem that generalizes the algorithm of Moerder and Calise [15], for discrete-time systems, is given (we call it the MC algorithm). The MC algorithm has a proof of convergence to a local optimum only, and it needs other algorithms for computing initial stabilizing SOF.

In Section 2 we formulate the problem and give some useful lemmas (without a proof). In Section 3, we introduce the randomized algorithm for the problem of

7.A comparison between the RS and the MC algorithms, as well as the performance of the hybrid algorithm, for real-life systems, is provided.

LQR via SOF for discrete-time LTI systems. Section 4 is devoted to the

deterministic algorithm for the problem. In Section 5, we give the results of the algorithms for some real-life systems. Finally, in Section 6 we conclude with

The reminder of the chapter is organized as follows:

equations (QMES) and thus can be applied to large systems.

3.Experience with the algorithm shows good quality of controllers (in terms of reduction of the LQR functional value with relatively small controller norms), 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

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. The Ray-Shooting Method was successfully applied recently also to the following hard complexity control problems for continuous-time

• Structured and structured-sparse SOFs (see [11])

(PID) controllers (see [13])

complexity.

some remarks.

**56**

for discrete-time systems.

for discrete-time systems.

• Robust control via SOF (see [14])

• LQR via SOF for continuous-time LTI systems (see [12])

systems:

*Control Theory in Engineering*

Let a discrete-time system be given by

$$\begin{cases} \boldsymbol{\omega}\_{k+1} = A \boldsymbol{\omega}\_k + B \boldsymbol{u}\_k, k = \mathbf{0}, \mathbf{1}, \dots \\ \boldsymbol{\chi}\_k = \mathbf{C} \boldsymbol{\omega}\_k \end{cases} \tag{1}$$

where *A* ∈ *<sup>p</sup>*�*p*, *B*∈ *<sup>p</sup>*�*<sup>q</sup>* , *C* ∈ *<sup>r</sup>*�*p*, and *x*<sup>0</sup> ∈ *<sup>p</sup>*. Let the LQR cost functional be defined by

$$J(\mathbf{x}\_0, \boldsymbol{\mu}) \coloneqq \sum\_{k=0}^{\infty} \left( \mathbf{x}\_k^T \mathbf{Q} \mathbf{x}\_k + \boldsymbol{\mu}\_k^T \mathbf{R} \boldsymbol{\mu}\_k \right), \tag{2}$$

where *Q* >0 and *R* >0. Let *uk* ¼ �*Kyk* be the SOF, and let *Ac*ℓð Þ *K* ≔ *A* � *BKC* denote the closed-loop matrix. Let denote the open unit disk, let 0 <*α* <1, and let *<sup>α</sup>* denote the set of all *z*∈ with j j *z* < 1 � *α* (where j j *z* is the absolute value of *z*). For a square matrix *Z*, we denote by *σ*ð Þ *Z* its spectrum. For any rectangular matrix *<sup>Z</sup>*, we denote by *<sup>Z</sup>*<sup>þ</sup> its Moore-Penrose pseudo-inverse. By k k*<sup>Z</sup> <sup>F</sup>* <sup>¼</sup> *trace ZTZ* � �<sup>1</sup> <sup>2</sup> we denote the Frobenius norm of *<sup>Z</sup>*, and by k k*<sup>Z</sup>* <sup>¼</sup> max *<sup>σ</sup> <sup>Z</sup>TZ* � � � � � � <sup>1</sup> <sup>2</sup> we denote the spectral norm of *Z*. By *LZ* and *RZ*, we denote the (left and right) orthogonal projections *I* � *Z*þ*Z* and *I* � *ZZ*<sup>þ</sup> on the spaces *Ker Z*ð Þ and *Ker Z*<sup>þ</sup> ð Þ, respectively. For a topological space X and a subset U ⊂ X, we denote by *int* ð Þ U the interior of U, i.e., the largest open set included in U. By U we denote the closure of U, i.e., the smallest closed set containing <sup>U</sup>, and by *<sup>∂</sup>*U ¼ U � *int* ð Þ <sup>U</sup> we denote the boundary of <sup>U</sup>. Let <sup>S</sup>*<sup>q</sup>*�*<sup>r</sup>* denote the set of all matrices *<sup>K</sup>* <sup>∈</sup> *<sup>q</sup>*�*<sup>r</sup>* such that *<sup>σ</sup>*ð Þ *Ac*<sup>ℓ</sup> <sup>⊂</sup> (i.e., stable in the discrete-time sense), and let <sup>S</sup>*<sup>q</sup>*�*<sup>r</sup>* denote the set of all matrices *<sup>K</sup>* <sup>∈</sup> *<sup>q</sup>*�*<sup>r</sup>* such that *σ*ð Þ *Ac*<sup>ℓ</sup> ⊂ *α*. If the last is nonempty, we say that *Ac*<sup>ℓ</sup> is *α*-stable and we call *α* the degree of stability. Let *<sup>K</sup>* <sup>∈</sup>S*<sup>q</sup>*�*<sup>r</sup> <sup>α</sup>* be given. Substitution of the SOF *uk* ¼ �*Kyk* ¼ �*KCxk* into (2) yields:

$$J(\mathbf{x}\_0, \mathbf{K}) = \sum\_{k=0}^{\infty} \mathbf{x}\_k^T \{ \mathbf{Q} + \mathbf{C}^T \mathbf{K}^T \mathbf{R} \mathbf{K} \mathbf{C} \} \mathbf{x}\_k. \tag{3}$$

Since *<sup>Q</sup>* <sup>þ</sup> *CTKTRKC* <sup>&</sup>gt;0 and *Ac*ℓð Þ *<sup>K</sup>* is stable, it follows that the Stein equation

$$P - A\_{c\ell}(K)^T P A\_{c\ell}(K) = Q + \mathbf{C}^T K^T R \mathbf{K} \mathbf{C} \tag{4}$$

has a unique solution *P* >0, given by

$$P(K) = \text{mat}\left(\left(I\_p \otimes I\_p - A\_{\epsilon\ell}(K)^T \otimes A\_{\epsilon\ell}(K)^T\right)^{-1} \cdot \text{vec}\left(Q + \mathbf{C}^T K^T \mathbf{R} \mathbf{K}\right)\right). \tag{5}$$

Substitution of (4) into (3) and using *xk* <sup>¼</sup> *Ac*ℓð Þ *<sup>K</sup> <sup>k</sup> x*<sup>0</sup> with the fact that *Ac*ℓð Þ *K* is stable leads to

$$\begin{split} J(\mathbf{x}\_{0},K) &= \sum\_{k=0}^{\infty} \mathbf{x}\_{k}^{T} \Big( P - \mathbf{A}\_{\varepsilon\ell}(K)^{T} P \mathbf{A}\_{\varepsilon\ell}(K) \Big) \mathbf{x}\_{k} \\ &= \sum\_{k=0}^{\infty} \mathbf{x}\_{0}^{T} \mathbf{A}\_{\varepsilon\ell}(K)^{Tk} \Big( P - \mathbf{A}\_{\varepsilon\ell}(K)^{T} P \mathbf{A}\_{\varepsilon\ell}(K) \Big) \mathbf{A}\_{\varepsilon\ell}(K)^{k} \mathbf{x}\_{0} \\ &= \mathbf{x}\_{0}^{T} P(K) \mathbf{x}\_{0} = \left\| P(K)^{\frac{1}{2}} \mathbf{x}\_{0} \right\|^{2}. \end{split}$$

Thus, when *<sup>x</sup>*<sup>0</sup> is known, we search for *<sup>K</sup>* <sup>∈</sup>S*q*�*<sup>r</sup> <sup>α</sup>* that minimizes the functional

$$J(\mathbf{x}\_0, K) = \mathfrak{x}\_0^T P(K) \mathfrak{x}\_0. \tag{6}$$

Lemma 2.2. We have:

*DOI: http://dx.doi.org/10.5772/intechopen.89319*

ð Þ *x*, *y x*∈*X*, *f x*ð Þ≤*y*≤ *y*<sup>0</sup>

Let *z*<sup>0</sup> ≔ *x*0, *y*<sup>0</sup>

Then

1.Let *A*, *B*,*X* be matrices with sizes *p* � *q*,*r* � *p*, *q* � *r*, respectively. Then

*<sup>∂</sup><sup>X</sup> trace AXB* ð Þ¼ *<sup>A</sup>TBT:*

2.Let *A*, *B*,*C*,*X* be matrices with sizes *p* � *q*,*r* � *r*, *q* � *p*,*r* � *q*, respectively.

*<sup>∂</sup><sup>X</sup> trace AXTBXC* <sup>¼</sup> *BTXATCT* <sup>þ</sup> *BXCA:*

The Ray-Shooting Method works as follows, for general function minimization: let *f x*ð Þ≥0 be a continuous function defined over some compact set *<sup>X</sup>* <sup>⊂</sup> *<sup>n</sup>*. Let ϵ >0 be given and assume that we want to compute *x*<sup>∗</sup> ∈*X* such that *y* <sup>∗</sup> ≔ *f x*ð Þ¼ <sup>∗</sup> min *<sup>x</sup>*∈*Xf x*ð Þ up to ϵ, i.e., to find ð Þ *x*, *y* in the set Sð Þ¼ ϵ fð Þ *x*, *y* j *x*∈*X*, *f x*ð Þ <sup>∗</sup> <sup>≤</sup> *<sup>y</sup>* <sup>¼</sup> *f x*ð Þ<sup>≤</sup> *f x*ð Þþ <sup>∗</sup> <sup>ϵ</sup>*:*g. Let *<sup>x</sup>*<sup>0</sup> <sup>∈</sup> <sup>X</sup> be given, let *<sup>y</sup>*<sup>0</sup> <sup>≔</sup> *f x*ð Þ<sup>0</sup> and let <sup>S</sup>ð Þ <sup>0</sup> <sup>¼</sup>

denote the search space, which is a subset of the epi-

between <sup>X</sup> and the level *<sup>y</sup>*0. Let <sup>L</sup>ð Þ <sup>0</sup> <sup>¼</sup> fð Þ *<sup>x</sup>*, *<sup>y</sup>* <sup>j</sup> *<sup>x</sup>*<sup>∈</sup> <sup>X</sup>, *<sup>y</sup>* <sup>¼</sup> <sup>0</sup> or *<sup>x</sup>*<sup>∈</sup> *<sup>∂</sup>*X, 0 <sup>≤</sup>*y*<sup>≤</sup> *<sup>y</sup>*0*:*g.

*tw*0, 0 ≤*t*≤ 1. We scan the ray and choose the largest 0 ≤*t*<sup>0</sup> ≤1 such that 1ð Þ � *t*<sup>0</sup> *z*<sup>0</sup> þ *<sup>t</sup>*0*w*<sup>0</sup> <sup>∈</sup>Sð Þ <sup>0</sup> (actually, we scan the ray from *<sup>t</sup>* <sup>¼</sup> 1 in equal-spaced points and take the first *t* for which this happens). We define *z*<sup>1</sup> ≔ ð Þ 1 � *t*<sup>0</sup> *z*<sup>0</sup> þ *t*0*w*<sup>0</sup> and update sets

<sup>L</sup>ð Þ<sup>1</sup> denote the updated sets. We continue the process similarly from *<sup>z</sup>*<sup>1</sup> <sup>∈</sup>Sð Þ<sup>1</sup> , and we define a sequence *zn* <sup>∈</sup>Sð Þ *<sup>n</sup>* , *<sup>n</sup>* <sup>¼</sup> 0, 1, …. Note that Sð Þ<sup>ϵ</sup> <sup>⊂</sup>Sð Þ *<sup>n</sup>*þ<sup>1</sup> <sup>⊂</sup> <sup>S</sup>ð Þ *<sup>n</sup>* for any *n* ¼ 0, 1, …, unless we have *zn* ∈Sð Þϵ for some *n* (in which the process is ceased).

Sð Þ*ε* . Note that shooting rays from the points of local minimum have positive probability to hit Sð Þ*ε* (under the following mild assumption), because any global minimum is visible from any local minimum. Moreover, for a given level of certainty, we hit Sð Þ*ε* in a finite number of iterations (see Remark 3.1 below). Practically, we may stop the algorithm if no improvement is detected within a window of 20% of the allowed number of iterations. The function need not be smooth or even continuous. It only needs to be well defined and measurable over the compact domain X, and Sð Þ*ε* should have non-negligible measure (i.e., should have some positive volume). Obviously, global minimum points belong to the boundary of the search space <sup>S</sup>ð Þ <sup>0</sup> , and actually such points are where the distance between the compact sets X � f g<sup>0</sup> and <sup>S</sup>ð Þ <sup>0</sup> in *<sup>n</sup>*þ<sup>1</sup> is accepted. This is essential for the efficiency of the Ray-Shooting Method, although we raised the search space

according to some distribution, and we define the ray as *z t*ð Þ ≔ ð Þ 1 � *t z*0þ

and note that *<sup>z</sup>*<sup>0</sup> <sup>∈</sup>Sð Þ <sup>0</sup> . Then, we choose *<sup>w</sup>*<sup>0</sup> in <sup>L</sup>ð Þ <sup>0</sup> randomly,

denote the cylinder enclosed

<sup>¼</sup> *<sup>z</sup>*1. Let <sup>S</sup>ð Þ<sup>1</sup> , <sup>D</sup>ð Þ<sup>1</sup> , and

*<sup>n</sup>*¼<sup>0</sup> converges (in probability) to a point in

*∂*

*Algorithms for LQR via Static Output Feedback for Discrete-Time LTI Systems*

**3. The randomized Ray-Shooting Method-based algorithm**

*∂*

graph of *<sup>f</sup>*. Let <sup>D</sup>ð Þ <sup>0</sup> <sup>¼</sup> ð Þ *<sup>x</sup>*, *<sup>y</sup> <sup>x</sup>*<sup>∈</sup> <sup>X</sup>, 0 <sup>≤</sup>*y*<sup>≤</sup> *<sup>y</sup>*<sup>0</sup>

One can show that the sequence f g *zn* <sup>∞</sup>

dimension from *n* to *n* þ 1.

**59**

<sup>S</sup>ð Þ <sup>0</sup> , <sup>D</sup>ð Þ <sup>0</sup> , and Lð Þ <sup>0</sup> by replacing *<sup>y</sup>*<sup>0</sup> with *<sup>y</sup>*1, where *<sup>x</sup>*1, *<sup>y</sup>*<sup>1</sup>

Let

$$\sigma\_{\max}(K) \coloneqq \max \left( \sigma(P(K)) \right). \tag{7}$$

Now, since *J x*ð Þ 0, *<sup>K</sup>* k k *<sup>x</sup>*<sup>0</sup> <sup>2</sup> <sup>≤</sup>*σmax*ð Þ *<sup>K</sup>* for any *<sup>x</sup>*<sup>0</sup> 6¼ 0, with equality in the worst case, therefore

$$\sup \sup\_{\mathbf{x}\_0 \neq 0} \left( \frac{\mathbf{J}(\mathbf{x}\_0, \mathbf{K})}{\left\| \mathbf{x}\_0 \right\|^2} \right) = \sup \mathbf{p}\_{\mathbf{x}\_0 \neq 0} \left( \frac{\left\| P(\mathbf{K})^{\frac{1}{2}} \mathbf{x}\_0 \right\|^2}{\left\| \mathbf{x}\_0 \right\|^2} \right) = \left\| P(\mathbf{K})^{\frac{1}{2}} \right\|^2 = \left\| P(\mathbf{K}) \right\| = \sigma\_{\max}(\mathbf{K}).$$

Thus, when *<sup>x</sup>*<sup>0</sup> is unknown, we search for *<sup>K</sup>* <sup>∈</sup>S*<sup>q</sup>*�*<sup>r</sup> <sup>α</sup>* , such that *σmax*ð Þ¼ *K* k k *P K*ð Þ is minimal. Note that if *λ* is an eigenvalue of *Ac*ℓð Þ *K* and *v* is a corresponding eigenvector, then (4) yields 1 � j j *<sup>λ</sup>* <sup>2</sup> <sup>¼</sup> *<sup>v</sup>* <sup>∗</sup> *<sup>Q</sup>*þ*CTKT* ð Þ *RKC <sup>v</sup> <sup>v</sup>* <sup>∗</sup> *P K*ð Þ*<sup>v</sup>* <sup>≥</sup> *<sup>v</sup>* <sup>∗</sup> *Qv <sup>v</sup>* <sup>∗</sup> *P K*ð Þ*<sup>v</sup>* <sup>&</sup>gt;0. Therefore, j j *<sup>λ</sup>* <sup>2</sup> <sup>≤</sup><sup>1</sup> � *<sup>v</sup>* <sup>∗</sup> *Qv <sup>v</sup>* <sup>∗</sup> *P K*ð Þ*<sup>v</sup>* <sup>&</sup>lt; 1, and thus, minimizing *<sup>σ</sup>max*ð Þ *<sup>K</sup>* results in eigenvalues that are getting closer to the boundary of . Since *α*, the degree of stability, is important to get satisfactory decay rate of the state to 0, and for disturbance rejection, we allow the user of the algorithms to determine *α*. Note that too high value of *α* might result in nonexistence of any SOF for the system or in complicating the search for a starting SOF. Higher values of *α* result in higher values of the optimal value of the LQR functional, i.e., higher energy consumption for decaying the disturbance *x*<sup>0</sup> to 0.

The functionals *J x*ð Þ 0,*K* and *σmax*ð Þ *K* are generally not convex since their domain of definition <sup>S</sup>*<sup>q</sup>*�*<sup>r</sup> <sup>α</sup>* (and therefore <sup>S</sup>*<sup>q</sup>*�*<sup>r</sup> <sup>α</sup>* ) is generally non-convex. Necessary conditions for optimality for continuous-time systems were given as three QMEs in [15–18]. Necessary and sufficient conditions for optimality for continuous-time systems, based on linear matrix inequalities (LMI), were given in [19–21]. 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.

In the sequel, we will use the following lemmas, given here without proofs. Lemma 2.1. We have:

1.The equation *AX* ¼ *B* has solutions if and only if *AA*þ*B* ¼ *B* or equivalently, if and only if *RAB* ¼ 0. In this case, the set of all solutions is given by

$$X = \mathcal{A}^+ \mathcal{B} + L\_A Z\_\*$$

where *Z* is arbitrary.

2.The equation *XA* ¼ *B* has solutions if and only if *BA*þ*A* ¼ *B* or equivalently, if and only if *BLA* ¼ 0. In this case, the set of all solutions is given by

$$X = BA^{+} + ZR\_{A},$$

where *Z* is arbitrary.

Lemma 2.2. We have:

Thus, when *<sup>x</sup>*<sup>0</sup> is known, we search for *<sup>K</sup>* <sup>∈</sup>S*q*�*<sup>r</sup>*

¼ *s upx*06¼<sup>0</sup>

Thus, when *<sup>x</sup>*<sup>0</sup> is unknown, we search for *<sup>K</sup>* <sup>∈</sup>S*<sup>q</sup>*�*<sup>r</sup>*

eigenvector, then (4) yields 1 � j j *<sup>λ</sup>* <sup>2</sup> <sup>¼</sup> *<sup>v</sup>* <sup>∗</sup> *<sup>Q</sup>*þ*CTKT* ð Þ *RKC <sup>v</sup>*

Let

therefore

*s upx*06¼<sup>0</sup>

j j *<sup>λ</sup>* <sup>2</sup> <sup>≤</sup><sup>1</sup> � *<sup>v</sup>* <sup>∗</sup> *Qv*

*x*<sup>0</sup> to 0.

domain of definition <sup>S</sup>*<sup>q</sup>*�*<sup>r</sup>*

Lemma 2.1. We have:

where *Z* is arbitrary.

where *Z* is arbitrary.

**58**

Now, since *J x*ð Þ 0, *<sup>K</sup>*

*Control Theory in Engineering*

*J x*ð Þ 0,*K* k k *x*<sup>0</sup> 2 ! *J x*ð Þ¼ 0, *<sup>K</sup> xT*

*P K*ð Þ<sup>1</sup> <sup>2</sup>*x*<sup>0</sup>

> k k *x*<sup>0</sup> 2

is minimal. Note that if *λ* is an eigenvalue of *Ac*ℓð Þ *K* and *v* is a corresponding

getting closer to the boundary of . Since *α*, the degree of stability, is important to get satisfactory decay rate of the state to 0, and for disturbance rejection, we allow the user of the algorithms to determine *α*. Note that too high value of *α* might result in nonexistence of any SOF for the system or in complicating the search for a starting SOF. Higher values of *α* result in higher values of the optimal value of the LQR functional, i.e., higher energy consumption for decaying the disturbance

The functionals *J x*ð Þ 0,*K* and *σmax*ð Þ *K* are generally not convex since their

conditions for optimality for continuous-time systems were given as three QMEs in [15–18]. Necessary and sufficient conditions for optimality for continuous-time systems, based on linear matrix inequalities (LMI), were given in [19–21].

However, algorithms based on these formulations are generally not guaranteed to converge, seemingly because of the non-convexity of the coupled matrix equations

In the sequel, we will use the following lemmas, given here without proofs.

and only if *RAB* ¼ 0. In this case, the set of all solutions is given by

and only if *BLA* ¼ 0. In this case, the set of all solutions is given by

1.The equation *AX* ¼ *B* has solutions if and only if *AA*þ*B* ¼ *B* or equivalently, if

*X* ¼ *A*þ*B* þ *LAZ*,

2.The equation *XA* ¼ *B* has solutions if and only if *BA*þ*A* ¼ *B* or equivalently, if

*X* ¼ *BA*<sup>þ</sup> þ *ZRA*,

*<sup>α</sup>* (and therefore <sup>S</sup>*<sup>q</sup>*�*<sup>r</sup>*

or inequalities, and when they converge, it is to a local optimum only.

� � �

0

B@

*<sup>α</sup>* that minimizes the functional

¼ k k *P K*ð Þ ¼ *σmax*ð Þ *K :*

*<sup>α</sup>* , such that *σmax*ð Þ¼ *K* k k *P K*ð Þ

*<sup>v</sup>* <sup>∗</sup> *P K*ð Þ*<sup>v</sup>* <sup>&</sup>gt;0. Therefore,

*<sup>α</sup>* ) is generally non-convex. Necessary

<sup>0</sup> *P K*ð Þ*x*0*:* (6)

*σmax*ð Þ *K* ≔ max ð Þ *σ*ð Þ *P K*ð Þ *:* (7)

2

� � � 2

k k *<sup>x</sup>*<sup>0</sup> <sup>2</sup> <sup>≤</sup>*σmax*ð Þ *<sup>K</sup>* for any *<sup>x</sup>*<sup>0</sup> 6¼ 0, with equality in the worst case,

1

CA <sup>¼</sup> *P K*ð Þ<sup>1</sup>

*<sup>v</sup>* <sup>∗</sup> *P K*ð Þ*<sup>v</sup>* <sup>≥</sup> *<sup>v</sup>* <sup>∗</sup> *Qv*

� � �

� � � 2

*<sup>v</sup>* <sup>∗</sup> *P K*ð Þ*<sup>v</sup>* <sup>&</sup>lt; 1, and thus, minimizing *<sup>σ</sup>max*ð Þ *<sup>K</sup>* results in eigenvalues that are

1.Let *A*, *B*,*X* be matrices with sizes *p* � *q*,*r* � *p*, *q* � *r*, respectively. Then

$$\frac{\partial}{\partial \mathbf{X}} \text{trace}(\mathbf{A}\mathbf{X}\mathbf{B}) = \mathbf{A}^T \mathbf{B}^T \mathbf{A}$$

2.Let *A*, *B*,*C*,*X* be matrices with sizes *p* � *q*,*r* � *r*, *q* � *p*,*r* � *q*, respectively. Then

$$\frac{\partial}{\partial \mathbf{X}} \text{trace} \left( \mathbf{A} \mathbf{X}^T \mathbf{B} \mathbf{X} \mathbf{C} \right) = \mathbf{B}^T \mathbf{X} \mathbf{A}^T \mathbf{C}^T + \mathbf{B} \mathbf{X} \mathbf{C} \mathbf{A}.$$
