Optimal control of discrete-time switched linear systems via continuous parameterization

TL;DR

This paper introduces a new optimization-based method for designing optimal controllers for discrete-time switched linear systems, reducing computational complexity by using auxiliary variables and non-smooth optimization techniques.

Contribution

It proposes a continuous parameterization approach that simplifies the switching law design, enabling more efficient computation compared to traditional methods.

Findings

Reduces computational complexity in switching law design

Uses auxiliary continuous input variables for optimization

Demonstrates effectiveness through simulations

Abstract

The paper presents a novel method for designing an optimal controller for discrete-time switched linear systems. The problem is formulated as one of computing the discrete mode sequence and the continuous input sequence that jointly minimize a quadratic performance index. State-of-art methods for solving such a control problem suffer in general from a high computational requirement due to the fact that an exponential number of switching sequences must be explored. The method of this paper addresses the challenge of the switching law design by introducing auxiliary continuous input variables and then solving a non-smooth block-sparsity inducing optimization problem.

Tables

Table 1. Table 1:

Threshold on ${\epsilon }_i$	${10}^{-5}$	${10}^{-7}$	${10}^{-8}$	${10}^{-10}$	$0$
$\%$ examples	$100\%$	$98\%$	$97\%$	$96\%$	$83\%$

Table 2. Table 2:

Threshold	${10}^{-2}$	${10}^{-5}$	${10}^{-7}$	${10}^{-8}$	${10}^{-10}$	$0$
$\%$ examples	$100\%$	$96\%$	$93\%$	$92\%$	$90\%$	$81\%$

Equations

$$x(k+1)=A_{\sigma (k)}x(k)+B_{\sigma (k)}u(k),x(0)=x_0$$

$$\begin{array}[]{l}J_{0}(x_{0},u(\cdot),\sigma(\cdot))=\dfrac{1}{2}x^{\textsf{T}}(N)\Psi x(N)\\ \hfill+\dfrac{1}{2}\sum\limits_{\begin{subarray}{c}k=0\end{subarray}}^{N-1}\big{[}x^{\textsf{T}}(k)Q(k)x(k)+u^{\textsf{T}}(k)R(k)u(k)\big{]}\end{array}$$

$$\displaystyle\min_{u(\cdot),\sigma(\cdot)}J_{0}\big{(}x_{0},u(\cdot),\sigma(\cdot)\big{)}$$

$$J_k^{\ast }(x)$$

$$\text{s.t. }x(t+1)=A_{\sigma (t)}x(t)+B_{\sigma (t)}u(t),x(k)=x$$

$$t=k,\dots ,N-1.$$

$$J_{k}^{*}(x)=\min_{\begin{subarray}{c}(u,i)\in\mathbb{R}^{n_{u}}\times\Omega\end{subarray}}\big{[}\ell(x,u,k)+J^{*}_{k+1}(A_{i}x+B_{i}u)\big{]}$$

$$J_k^{\ast }(x)=x^{\text{T}}{P}^{\ast }(k)x\forall x\in {\mathbb{R}}^n,\forall k\in \left\{0,\dots ,N\right\},$$

$$P^{\ast }(k)={\rho }_{{\sigma }^{\ast }(k),k}(P^{\ast }(k+1)),P^{\ast }(N)=\Psi$$

$${\rho }_{i,k}(P)=$$

$$\displaystyle-A_{i}^{\textsf{T}}PB_{i}\big{(}R(k)+B_{i}^{\textsf{T}}PB_{i}\big{)}^{-1}B_{i}^{\textsf{T}}PA_{i}.$$

$$u^{*}(k)=-K_{\sigma^{*}(k),k}\big{(}P^{*}(k+1)\big{)}x^{*}(k)$$

$$\sigma^{*}(k)\in\underset{i\in\Omega}{\text{argmin}}\big{[}x^{*}(k)^{\textsf{T}}\rho_{i,k}(P^{*}(k+1))x^{*}(k)\big{]}.$$

$$J_k^{\ast }(x)=\underset{P\in {\mathcal{H}}_k}{\min }{x}^{\text{T}}Px$$

$$\left\{\begin{aligned} &\mathcal{H}_{N}=\{\Psi\}\quad\\ &\mathcal{H}_{k}=\big{\{}\rho_{i,k}(P):P\in\mathcal{H}_{k+1},i\in\Omega\big{\}},\>k=N-1,\ldots,0.\end{aligned}\right.$$

$${\sigma }^{\ast }(k)\in \underset{i\in \Omega ,P\in {\mathcal{H}}_{k+1}}{\mathrm{argmin}}{x}^{\ast }(k{)}^{\top }{\rho }_{i,k}(P)x^{\ast }(k)$$

$$x(k+1)=A_{i}x(k)+B_{i}u(k)+f_i(k),\forall i\in \Omega$$

$$\begin{array}[]{l}\mathcal{S}=\big{\{}v=(v_{1}^{\textsf{T}},\cdots,v_{q}^{\textsf{T}})^{\textsf{T}}\in\mathbb{R}^{qn},\mbox{ with }v_{i}\in\mathbb{R}^{n},i\in\Omega:\\ \hfill\exists j\in\Omega:v_{j}=\textbf{0}_{n,1}\big{\}}.\end{array}$$

$$\displaystyle\mathcal{S}=\Big{\{}v=(v_{1}^{\textsf{T}},\cdots,v_{q}^{\textsf{T}})^{\textsf{T}}\in\mathbb{R}^{qn},\mbox{ with }$$

$$\displaystyle\Big{.}\prod_{i=1}^{q}\left\|v_{i}\right\|=0\Big{\}}$$

$$V(x_0,$$

$$\left[\begin{array}[]{c}x(k+1)\\ \vdots\\ x(k+1)\end{array}\right]=\left[\begin{array}[]{c}A_{1}\\ \vdots\\ A_{q}\end{array}\right]x(k)+\left[\begin{array}[]{c}B_{1}\\ \vdots\\ B_{q}\end{array}\right]u(k)+\bar{f}(k)$$

$$\bar{A}={\left[\begin{array}{ccc}{A}_1^{\text{T}}& \cdots & A_q^{\text{T}}\end{array}\right]}^{\text{T}},\bar{B}={\left[\begin{array}{ccc}{B}_1^{\text{T}}& \cdots & B_q^{\text{T}}\end{array}\right]}^{\text{T}},$$

$$L_{n,q}={\left[\begin{array}{cc}{I}_n& \cdots I_n\end{array}\right]}^{\text{T}}={\text{1}}_q\otimes I_n,$$

$$\underset{u(\cdot ),\bar{f}(\cdot )}{\min }$$

$$\text{s.t. }{L}_{n,q}x(k+1)=\bar{A}x(k)+\bar{B}u(k)+\bar{f}(k)$$

$$\bar{f}(k)\in \mathcal{S},k=0,\dots ,N-1$$

$$x(0)=x_0.$$

$$\underset{u(\cdot )}{\min }$$

$$L_{n,q}x(k+1)=\bar{A}x(k)+\bar{B}u(k)+\bar{f}(k),$$

$$k=0,\dots ,N-1$$

$$x(0)=x_0,$$

$$\begin{array}[]{l}\bar{V}=V\big{(}x_{0},u(\cdot),\bar{f}(\cdot)\big{)}+\\ \sum_{k=0}^{N-1}\bar{\lambda}^{\textsf{T}}(k+1)\big{(}\bar{A}x(k)+\bar{B}u(k)+\bar{f}(k)-L_{n,q}x(k+1)\big{)}\end{array}$$

$$H_k$$

$$=\ell (x(k),u(k),k)+{\lambda }^{\text{T}}(k+1)\left(\bar{A}x(k)+\bar{B}u(k)+\bar{f}(k)\right)$$

$$\bar{V}=\frac{1}{2}x(N{)}^{\text{T}}$$

Full text

1 Introduction

Switched systems constitute a class of dynamic systems consisting of a finite number of subsystems which are activated one after another over time by a switching signal (12), (8), (7). In many cases, the switching signal is an external input which, together with the continuous input vector, can be used to control the behavior of the switched system. Examples of real processes that can be represented as switched systems with external discrete and/or continuous inputs are: autonomous vehicles, chemical processes, electrical circuits, etc. (see e.g., (12) and (8) for a background). The problem of controlling switched systems of this type therefore involves designing a continuous control law along with a switching sequence to achieve some performance specifications. We will be interested more specifically in optimal control design for discrete-time switched linear systems (SLS).

The existing approaches to the problem of optimal control for switched systems can be classified into two groups: the ones addressing the continuous-time case (13; 5; 9; 10; 11) and those pertaining to the discrete-time case (14), (6). Other classifications can be done with respect to the length of the control horizon (finite or infinite), the simplifying assumptions posed a priori on the structure of the discrete or the continuous input or the conceptual nature of the methods. One observation that can be made is the following. Considering the quadratic optimal control of SLS in continuous-time where no structure is imposed in advance on the discrete and continuous control policies, there is so far no exact solution for both finite and infinite control horizons. In discrete-time an exact solution has been derived in (14) for the case of finite control horizon. However direct numerical implementation proves to be so expensive that it is not affordable in practice. Therefore a relevant question of major importance is how to develop some suboptimal strategies that would be much less expensive while still being close to optimality. Some relaxations have been discussed in (14; 6) for this purpose. However the resulting control algorithms still suffer from an exponential demand in storage capacity.

In this paper, we focus on the quadratic optimal control of switched linear systems. The discussion here is restricted to finite-horizon problems but the envisioned ultimate goal is to extend it to infinite horizon. We first observe that the origin of the huge complexity associated with computing the solution to this optimization problem is the presence of discrete variables. Therefore a key idea in our approach is to parameterize the discrete inputs by auxiliary control variables which are continuous. With this reformulation, the SLS optimal control problem becomes completely continuous but nonconvex. We then discuss some convex relaxation strategies. More precisely, we propose solving a sequence of convex problems in order to estimate the auxiliary control variables. The whole process yields an implementation which is shown in simulation to coincide statistically very often with the true optimal control. The advantage of the proposed computational scheme is that it has only a polynomial complexity. Moreover it requires the same storing capacity (for the sequence of positive semidefinite matrices generated by the Riccati recursion) as the solution to the simpler linear quadratic optimal problem.

The structure of this paper is as follows. In Section 2, we formulate the switched quadratic control problem. We discuss the expression of the optimal solution and the associated complexity issue. In Section 3 we present the new continuous parameterization of the discrete control variable and develop a four-steps algorithm for solving it approximately. A numerical illustration is provided in Section 4. Finally, some concluding remarks are given in section 5.

2 Problem Formulation

2.1 Switched linear systems

We consider a discrete-time switched linear system (SLS) described by

[TABLE]

where $k\in\mathbb{Z}_{+}$ is a time index, $x(k)\in\mathbb{R}^{n}$ and $u(k)\in\mathbb{R}^{n_{u}}$ are respectively the state and the continuous input at time $k$, $x_{0}$ is the initial state; $\sigma(k)\in\Omega\triangleq\left\{1,\ldots,q\right\}$ refers to the value taken by the switching signal (also called here the discrete input) at time $k$. $\Omega$ is a finite set collecting the indices of the different subsystems of the SLS (1). For any $i\in\Omega$, the pair $(A_{i},B_{i})\in\mathbb{R}^{n\times n}\times\mathbb{R}^{n\times n_{u}}$ of matrices is associated with the subsystem $i$. Throughout the paper, we use notations of the type $u(\cdot)$ to designate the entire sequence $\left\{u(k)\right\}$.

It is important to note that the switching signal $\sigma(\cdot)\triangleq\left\{\sigma(k)\right\}$ is viewed here as an external input. For simplicity, we assume that $\sigma(k)$ can be selected freely in $\Omega$ without any constraint. The control problem of interest is that of designing jointly a continuous control sequence $u(\cdot)$ and a discrete control sequence $\sigma(\cdot)$ so as to minimize a certain performance index over a finite time horizon. For this purpose it will be assumed throughout the paper that the SLS (1) is stabilizable.

2.2 Switched optimal quadratic problem

We consider a quadratic performance index associated with system (1) in the form

[TABLE]

where $N$ denotes the control horizon, $\left\{R(k)\right\}_{k=0}^{N-1}\subset\mathbb{R}^{n_{u}\times n_{u}}$ is a given sequence of symmetric positive definite matrices; $\left\{Q(k)\right\}_{k=0}^{N}\subset\mathbb{R}^{n\times n}$ with $Q(N)=\Psi$ represents a sequence of positive semidefinite matrices. The problem of interest in this paper is stated as follows.

Problem 1

Given the performance matrices $\left\{Q(k)\right\}_{k=0}^{N}$ and $\left\{R(k)\right\}_{k=0}^{N-1}$, find a continuous input sequence $u(\cdot)$ and a discrete input sequence $\sigma(\cdot)$ that minimize the performance index (2) subject to the switched system equation (1). In more formal terms, this is equivalent to solving the optimization problem

[TABLE]

We start by observing that the solution to problem (3) can be mathematically characterized in a quite straightforward way using, for example, the Bellman optimality principle (4). However a direct implementation of the optimal control suffers from an exponential complexity in both storing and computational resources. The goal of this paper is to discuss some alternative formulations and corresponding solutions of the quadratic optimal control problem for SLS so as to yield more efficient implementations.

Let us first characterize the solution to problem (3). For this purpose, let $J_{k}(x,u(\cdot),\sigma(\cdot))$ denote the performance index corresponding to the situation when the system starts in an arbitrary state $x\in\mathbb{R}^{n}$ at time $k$ and evolves under the action of the inputs $u(\cdot)$ and $\sigma(\cdot)$. Introduce the function $J_{k}^{*}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ defined by

[TABLE]

$J_{k}^{*}(x)$ is called the value function or the cost-to-go. It is interesting to note from (2) that $J_{N}^{*}(x)=x^{\textsf{T}}\Psi x$ for all $x\in\mathbb{R}^{n}$. Equipped with the notation (4), the Bellman principle of optimality (4) can be expressed as

[TABLE]

where $\ell(x,u,k)=\dfrac{1}{2}\left[x^{\textsf{T}}Q(k)x+u^{\textsf{T}}R(k)u\right]$ is the running cost.

Theorem 1

Consider the control problem (3) and let the functions $J_{k}^{*}$ be defined as in (4). Denote with $u^{*}(\cdot)$, $\sigma^{*}(\cdot)$, $x^{*}(\cdot)$ respectively the optimal continuous input, discrete input and continuous state. Then the following statements hold.

The value function $J_{k}^{*}$ is quadratic and can be written as

[TABLE]

where $\left\{P^{*}(k)\right\}$ is a sequence of matrices generated recursively backward in time according to the Ricatti recursion defined by

[TABLE]

with

[TABLE] 2. 2.

The optimal continuous input is given by

[TABLE]

where K_{i,k}(P)=\big{(}R(k)+B_{i}^{\textsf{T}}PB_{i}\big{)}^{-1}B_{i}^{\textsf{T}}PA_{i}. 3. 3.

The optimal discrete input $\sigma^{*}(\cdot)$ is given by

[TABLE]

{pf}

The proof follows by a simple backward induction exploiting the Bellman optimality equation (5). It is therefore omitted.

2.3 Implementation of the optimal control

In order to implement the control law (9)-(10), we need to fully compute offline the sequence $\{P^{*}(k)\}_{k=0}^{N}$ by the Ricatti recursion (7). Note that this recursion depends on the optimal switching signal $\sigma^{*}(\cdot)$ which in turn depends on the optimal continuous state $x^{*}(\cdot)$ as can be seen from Eq. (10). Unfortunately the optimal state is not available offline. This is a source of a major difficulty. To cope with this challenge an elegant solution has been developed in (14). The authors first showed that the value function given in (6) can be re-expressed as

[TABLE]

where $\left\{\mathcal{H}_{k}\right\}$ is a sequence of sets of symmetric positive semidefinite matrices defined by111Note that in the setting of (14) the weighting matrices $Q$ and $R$ depend on the subsystem index, not on time. As a consequence $\rho$ is indexed there only by $i$.

[TABLE]

The sets $\left\{\mathcal{H}_{k}\right\}$ collect indeed all the possible values of the Ricatti sequence (7) for all admissible switching signals. Since the definition of the sets $\left\{\mathcal{H}_{k}\right\}$ is now freed from the dependence on the continuous state, they can be computed offline and stored. Once this is done the optimal control law can be obtained online by applying

[TABLE]

along with (9).

It turns out that the trick of (14) makes the optimal control (9)-(10) implementable. This is done however at the price of a huge complexity. In effect, the cardinality of the sets $\mathcal{H}_{k}$ grows exponentially fast with respect to the control horizon $N$. For example, the cardinality of $\mathcal{H}_{k}$ is about $q^{N-k}$. Indeed the above implementation requires an exponential load in terms of both computational and storage resources. This complexity affects all the steps of the implementation: offline computation and storage of the sets $\left\{\mathcal{H}_{k}\right\}$ and online reading and search for the optimal discrete control by Eq. (12) over the sets $\left\{\mathcal{H}_{k}\right\}$. To reduce the complexity, a suboptimal solution is discussed in (14). But it is fair to observe that some shortcomings still persist. First, the proposed procedure does not alleviate the off-line computational load ; it only allows for a saving of the necessary storage capacity. Moreover the cardinality reduction algorithm has still to test all the elements in the sets $\mathcal{H}_{k}$ hence resulting in an exponential complexity.

This complexity restricts the applicability of the solution of (14) to the control of switched systems with small number of subsystems and small control horizons. Noting that exponential complexity is generated by the presence of discrete variables in the optimization problem, we discuss here a new approach which relies on a continuous parameterization of the switching sequence.

3 Proposed solution

The proposed design method relies on two main ideas:

• •

A continuous parameterization of the discrete control variable $\sigma(\cdot)$. This consists in replacing the discrete input in the SLS (1) with continuous variables called auxiliary control variables. Consequently, problem (3) can be written as a constrained optimization problem in only continuous variables hence getting rid of its combinatorial feature.

• •

A nonsmooth block-sparsity-inducing optimization involving the auxiliary control variables. The purpose of this is to enforce their expected structure as will be described next.

3.1 Continuous parameterization of the discrete input

A starting point is to notice that the SLS dynamics in (1) can be written in the form

[TABLE]

where $f_{i}(k)=(A_{\sigma(k)}-A_{i})x(k)+(B_{\sigma(k)}-B_{i})u(k)$. Each $f_{i}(k)$ can be interpreted as the difference between the state of the SLS at time $k+1$ under $\sigma(k)$ and the state that would have been obtained if $\sigma(k)$ was set equal to $i$. It follows that for any time instant $k$, $f_{\sigma(k)}(k)=0$.

From now on let us forget about the explicit expressions of the $f_{i}(k)$’s and view them just as unknown control variables satisfying the following constraint: for all $k$, there exists $j=j(k)\in\Omega$ such that $f_{j}(k)=0_{n,1}$ where $0_{n,1}$ denotes a $n$-dimensional vector with all entries equal to zero. By letting $\bar{f}(k)=[\begin{matrix}f_{1}(k)^{\textsf{T}}&\cdots&f_{q}(k)^{\textsf{T}}\end{matrix}]^{\textsf{T}}\in\mathbb{R}^{qn}$, the above constraint can be written as $\bar{f}(k)\in\mathcal{S}$ with $\mathcal{S}$ being a subset of $\mathbb{R}^{qn}$ defined by

[TABLE]

There are many other equivalent ways of representing the set $\mathcal{S}$. One of those which are smooth is the following

[TABLE]

for any vector norm $\left\|\cdot\right\|$ on $\mathbb{R}^{n}$. The variable $\bar{f}(k)$ is then regarded as an auxiliary control variable to be computed along with the continuous input variable. Based on this parameterization we now restate the control problem as follows.

Problem 2

Given the matrices $\left\{Q(k)\right\}_{k=0}^{N}$ and $\left\{R(k)\right\}_{k=0}^{N-1}$, find a continuous input sequence $u(\cdot)$ and an auxiliary input sequence $\bar{f}(\cdot)$ that minimize the performance index

[TABLE]

subject to system (13) written in expanded form as

[TABLE]

and the constraint $\bar{f}(k)\in\mathcal{S}$.

To write this in a more compact form, define

[TABLE]

where $\textbf{1}_{q}=[\begin{matrix}1&\cdots&1\end{matrix}]^{\textsf{T}}$ is a vector of $q$ ones, $I_{n}$ is the $n$-dimensional identity matrix and $\otimes$ denotes the Kronecker product. Then problem (2) can be written as a constrained optimization problem as follows

[TABLE]

It is important to note that problems (3) and (17) are equivalent. But contrary to (3), problem (17) is an optimization problem in which all the decision variables are continuous. One can therefore hope for designing an algorithm that solves it in polynomial time. A remaining challenge to deal with is the non convexity of (17).

3.2 An algorithm in four stages

We now ask the question of how to compute numerically the solution to problem (17). Because of the non-convex constraint $\bar{f}(k)\in\mathcal{S}$, the problem is not convex. For efficiency of solving we therefore need to find a convex relaxation. We will discuss a block-sparsity inducing optimization technique for that. The global computational procedure can be decomposed into four stages which are described next.

3.2.1 Stage (a)

In this first stage, the sequence $\left\{\bar{f}(k)\right\}_{k=0}^{N-1}$ is considered as a (unknown) parameter. Hence the criterion $V$ is minimized with respect to $u(\cdot)$ only, i.e., we solve

[TABLE]

where the constraint on $\left\{\bar{f}(k)\right\}$ has been removed. This is a convex program. One difference with the classical linear quadratic control problem is that the dynamic matrix $\bar{A}$ in (18) is rectangular rather than square. Also the state $x(k+1)$ is repeated here $q$ times on left hand side of the dynamics equation. Denote with $V_{1}(x_{0},\bar{f}(\cdot))$ the resulting optimal functional, i.e., V_{1}(x_{0},\bar{f}(\cdot))=\min_{u(\cdot)}V\big{(}x_{0},u(\cdot),\bar{f}(\cdot)\big{)}.

Using the method of Lagrange multipliers, construct the Lagrangian $\bar{V}$ by embedding the constraints of (18) in the cost functional. This yields

[TABLE]

where $\left\{\bar{\lambda}(k)\right\}$, with $\bar{\lambda}(k)=[\begin{matrix}\lambda_{1}(k)^{\textsf{T}}&\cdots&\lambda_{q}(k)^{\textsf{T}}\end{matrix}]^{\textsf{T}}$, is a sequence of Lagrange multipliers. Introduce now the discrete time Hamiltonian $H$ associated with the system (16) and the performance index (15) defined by

[TABLE]

With this short-hand notation the extended cost (19) takes the form

[TABLE]

Considering the minimization of $\bar{V}$ with respect to $u$, let us look at the effect of an elementary variation $du$ of $u$ and of the induced change $dx$ in the state. These together induce a variation $d\bar{V}$ of $\bar{V}$, expressible by

[TABLE]

Note that $dx(0)=0$ since the initial state is fixed. At the optimal value of $\bar{V}$ we must have $d\bar{V}=0$ as a consequence of any variation $du$ of the input and any subsequent variation $dx$ of the state. This implies that all terms in the previous equation must be set to zero. Thus for all $k=1,\ldots,N-1$, $\bar{\lambda}(k)$ must satisfy

[TABLE]

Similarly, we must impose

[TABLE]

It follows that the optimal continuous input can be expressed as

[TABLE]

We can rewrite (16) and (23) respectively as follows

[TABLE]

with initial condition $x(0)=x_{0}$ and

[TABLE]

with the final constraint $L_{n,q}^{\textsf{T}}\bar{\lambda}(N)=\Psi x(N)$.

Eqs (25) and (26) form a system of linear equations that characterizes completely the solution to (18).

3.2.2 Stage (b)

Now, we seek the minimal value of $V_{1}(x_{0},\bar{f}(\cdot))$ with respect to $\bar{f}(\cdot)$. Since for any $k$ the solution $\bar{f}(k)$ is expected to live in the set $\mathcal{S}$ defined in (14), the optimization problem is as follows

[TABLE]

As mentioned earlier, this problem is not convex. In order to keep polynomial complexity, we need to find a convex relaxation of it.

To do so, let us denote with $\chi_{\mathcal{S}}:\mathbb{R}^{nq}\rightarrow\mathbb{R}\cup\left\{+\infty\right\}$ the indicator function of $\mathcal{S}$ defined from $\mathbb{R}^{nq}$ to real extended line by $\chi_{\mathcal{S}}(z)=0$ if $z\in\mathcal{S}$ and $\chi_{\mathcal{S}}(z)=+\infty$ otherwise. Then (27) is equivalent

[TABLE]

To find a convex relaxation of the terms $\chi_{\mathcal{S}}(\bar{f}(k))$, we view each vector $\bar{f}(k)$ as being relatively block-sparse in the sense that it must admit at least one subvector $f_{j}(k)$ which is equal to zero. With this in mind we replace $\chi_{\mathcal{S}}(\bar{f}(k))$ with a nonsmooth convex function, $\sum_{i=1}^{q}w_{i}(k)\left\|f_{i}(k)\right\|_{2}$ where the $w_{i}(k)$’s denote some positive weights. As is suggested by a certain number of results (see e.g., (3; 2)), minimizing such a function enjoys the nice property that it is able to promote block sparsity hence yielding a vector $\bar{f}(k)$ with some of its subvectors $f_{i}(k)$ potentially equal or close to zero. The role of the weights $\left\{w_{i}(k)\right\}$ is to discriminate between the different subvectors of $\bar{f}(k)$. They can be selected, for example, iteratively as follows: solve the problem with all weights set to one and based on the resulting solution, retune the weights through a simple rule of the form $w_{i}(k)=1/(\|\hat{f}_{i}(k)\|_{2}+\epsilon)$ for some small number $\epsilon>0$. For better numerical stability, one can consider normalizing a posteriori the weights along the $i$ dimension.

Finally we formulate the following convex optimization problem:

[TABLE]

For writing simplicity, we have not given here the explicit expression of $V_{1}(x_{0},f(\cdot))$. The derivation of such an expression follows from straightforward algebraic calculations departing from the system of linear equations (25)-(26).

3.2.3 Stage (c)

Let \big{\{}\hat{\bar{f}}(k)\big{\}} with $\hat{\bar{f}}(k)=[\begin{matrix}\hat{f}_{1}^{\top}(k)&\cdots&\hat{f}_{q}^{\top}(k)\end{matrix}]^{\top}$ denote a solution to problem (28). Since this is obtained under some relaxation of the initial problem, there is no guarantee that $\hat{\bar{f}}(k)$ will lie in the set $\mathcal{S}$ for all $k$. Hence, we need to project those $\hat{\bar{f}}(k)$ onto $\mathcal{S}$ in order to determine the switching sequence. The projection we used consists in selecting the discrete input to be the index $i\in\Omega$ such that $\hat{f}_{i}$ has the minimum norm among all, i.e.,

[TABLE]

with $\left\|\cdot\right\|$ denoting the vector $2$-norm. In some sense, this corresponds to forcing to zero the $\hat{f}_{i}$ having the minimum norm.

3.2.4 Stage (d)

Once we have computed the switching sequence $\{\hat{\sigma}(k)\}_{k=0}^{N-1}$ offline, the SLS optimal control problem reduces to that of a linear time-varying system with matrices defined by $(A_{\hat{\sigma}(k)},B_{\hat{\sigma}(k)})$. The solution of such a problem with respect to the quadratic cost (2) can be determined as in the conventional case (1). To obtain it, we can just apply (7)-(8) with $\sigma^{*}(\cdot)$ replaced by $\hat{\sigma}(\cdot)$ to generate offline a sequence of matrices $\{\hat{P}(k)\}$. By storing this single Ricatti sequence, the final continuous and discrete inputs are then, similarly as in (9)-(10), selected online as

[TABLE]

with $\hat{x}(\cdot)$ denoting the associated state trajectory. This means that the discrete input is recalculated online.

As it turns out, the implementation proposed in this paper has a polynomial complexity. Moreover it requires storing only a single sequence of matrices $\{\hat{P}(k)\}$ just as in the solution of the linear quadratic problem (single linear subsystem). In comparison with (14), the gain on the memory demand is significant. However in its current version, the proposed implementation is not guaranteed to yield the optimal control. In view of the applied control policy (30)-(31), the sole objective of the procedure described above for solving Problem 2 is the computation of the Ricatti sequence $\{\hat{P}(k)\}$.

4 Results

Note that for all the following tests, the weights $\gamma_{1}$ and $\gamma_{2}$ are taken equals to $1$.

4.1 A statistical test

In this section we challenge the capability of the proposed method to obtain the solution to the optimal quadratic control problem for switched systems of the form (1). For this purpose $100$ examples of two-dimensional switched systems, with each composed of $q=2$ subsystems, are generated at random using the MATLAB function drss. The initial state is also sampled from a Gaussian distribution $\mathcal{N}(0,20I_{2})$. For each of these examples, the optimal control problem is that of finding the discrete and continuous inputs to minimize a finite horizon performance index of the form (2) with $N=15$ and with $Q(k)=I_{2}$ and $R(k)=1$ for all $k$.

For the sake of comparison, we also implement the optimal control law as described in (14) and recalled above in Eqs (11)-(12). This is possible here since the number of subsystems and the control horizon are small. As a matter of fact the cardinality of $\mathcal{H}_{0}$ in (11) for the current experiment is about $2^{15}=32768$, which is computationally affordable on a standard computer. Let $J_{i}^{\text{opt}}$ and $J_{i}^{o}$, $i=1,\ldots,100$, denote respectively the optimal index and the value yielded by our method. For each of the $100$ examples, define the following relative error as an empirical measure of the distance to optimality

[TABLE]

The table below displays the distribution of $\varepsilon_{i}$ in terms of probabilities of the type $\operatorname*{Prob}(\varepsilon_{i}\leq\alpha)$ with $\alpha$ being a threshold taking values in $\left\{10^{-10},10^{-8},10^{-7},10^{-5}\right\}$. It turns out that the solution achieved by the proposed method either coincides with the optimal index or lies generally in a small neighborhood of it. This suggests that for generic systems, the implemented strategy discussed in this paper has the potential of recovering the optimal control and this, at a much affordable price.

Repeating a similar experiment as above with this time examples of switched systems with state dimension $n=3$ and number of modes $q=3$ yield the results reported in Table 2. Note that for this last experiment the control horizon has been reduced to $N=10$ in order to alleviate the computational load associated with the computation of the exact optimal solution. A little degradation of the results can be observed in Table 2 in comparison with the results given in Table 1. This may be due to numerical artefacts as a result of increased number of decision variables. The approximate performance index is still very close to the optimal one.

4.2 Illustration of performance on a single example

For illustration purpose, let us now focus on a single switched system. The considered example is in the form (1) and consists of two linear subsystems with matrices defined by

[TABLE]

It can be observed that none of the individual subsystems is stable. In this experiment, the control horizon is set to $N=15$; the initial state is $x_{0}=[\begin{matrix}1&2\end{matrix}]^{\top}$ and the weighting matrices of the performance index are defined as in Section 4.1. Applying the proposed method on this example yields the results presented in Figure 1. It turns out that the obtained discrete input is equal to the optimal one except at the two time instants $t=10$ and $t=14$. However the impact of this difference is negligible on the performance index $J_{0}$ since we still get a relative error of $4.03\times 10^{-9}$. This is because the errors occur at a time when the state has almost already converged to zero as shown by Figure 1-(b). We can even conjecture that small amplitude of the state makes it difficult to recover the optimal discrete input by (31).

5 Conclusion

In this paper, we studied the discrete-time quadratic optimal control problem for switched systems on a finite time horizon. Based on a continuous parameterization of the discrete input, we proposed an approach that is able to yield the optimal solution in polynomial time with respect to the length of control horizon. Moreover the proposed algorithm appears to be cheaper than most existing methods in terms of computational load and storage resources. Future research will focus on analyzing the properties of this method.