Computation of robust positively invariant set based on direct data-driven approach

Qian Du; Hongli Yang

doi:10.20517/ces.2024.76

Download PDF

Research Article | Open Access | 29 Nov 2024

Computation of robust positively invariant set based on direct data-driven approach

Views: 72 | Downloads: 15 | Cited:

0

Qian Du

,

Hongli Yang

Complex Eng Syst 2024;4:24.

10.20517/ces.2024.76 | © The Author(s) 2024.

Author Information

Article Notes

Cite This Article

Abstract

A direct data-driven approach for computing the robust positively invariant sets of a discrete-time system is presented in this study. By transforming the invariance conditions into a set of tractable linear matrix inequality problems and combining them with a semidefinite programming problem, we maximize the volume of invariant sets without violating state constraints. Based on two equivalence conditions of invariance, we investigate two algorithms using the one-step method to maximize the volume of the invariant sets. Ultimately, we opted for Algorithm 1, which is more succinct and effective. To further reduce conservatism, we propose an iterative algorithm based on Algorithm 1. The effectiveness of the proposed algorithm is verified through numerical examples.

Keywords

Direct data-driven, robust positively invariant set, linear discrete-time system, semidefinite programming

Download PDF 0 2

1. INTRODUCTION

A set is defined as a robust positively invariant (RPI) set, if every initial state contained within it, the trajectory of the system's state remains confined to that set, regardless of any disturbances or uncertainties in its parameters ^[1]. In the domain of automatic control, the significance of RPI sets is self-evident. The theories and computational approaches of RPI sets have not only attracted significant interest in the academic community but also demonstrated extensive application values in multiple practical fields such as system stability analysis, controller design, and nonlinear system theory ^[1-6]. Moreover, RPI sets offer us a powerful tool for evaluating and resolving a series of issues closely related to external unknown disturbances.

For linear discrete systems, several approaches have been proposed to compute the RPI sets ^[7-10]. These are model-based approaches that assume the system is usable. However, achieving precise and dependable models in practical applications is challenging, and imprecise models can result in a loss of system invariance and breach associated constraints ^{[11, 12]}, potentially leading to unstable system operation or suboptimal controller performance. To address these limitations, direct data-driven approaches have emerged as a viable alternative, eliminating the need for model identification ^[12]. As we all know, the two important forms of invariant sets are polyhedral and ellipsoidal sets ^[4]. Ellipsoidal sets may be more conservative than polyhedral sets. Since polyhedral sets offer more flexible and complex representations, they have an advantage over ellipsoidal sets both in theory and practice ^[1]. Additionally, polyhedral representations naturally capture physical constraints on state and control variables ^[3]. Due to these factors, this paper emphasizes the study of polyhedral sets.

In this paper, we propose a direct data-driven method ^{[11, 13]} for computing polyhedral RPI sets. We assume that the RPI set is a 0-symmetric convex polyhedron of predetermined complexity. The method presented in ^[14] involves a relatively high number of algorithm optimization variables and linearized matrix inequality constraints, resulting in slow running speeds. Our method expands upon and optimizes the approach recently developed by Mejari $$ et\; al $$^[14]. We have derived a new invariance condition and developed an iterative algorithm based on this condition. In addressing the nonlinear problem of the invariance condition, we ingeniously utilize relevant variable transformations and appropriate scaling techniques. Subsequently, a series of LMIs ^[15] are employed to precisely represent the constraints and invariance requirements of the system. Finally, we resort to the semidefinite programming (SDP) problem for solution, aiming to maximize the volume of the RPI sets. Notably, our method does not require an accurate system model or system identification steps; it merely requires a state trajectory composed of a finite number of data samples as input. Furthermore, the number of optimization variables and LMIs in this iterative algorithm is less than that in ^[14], thereby reducing running time and operating costs to some extent.

The primary contribution of this paper lies in the introduction of a novel data-driven methodology tailored for computing systems equipped with RPI sets. This approach holds significant relevance, particularly within the realms of control system design and analysis. A standout feature of our method is its reliance on a limited dataset, bypassing the need for comprehensive mathematical model knowledge-a pivotal advantage in numerous practical applications where precise mathematical models are elusive or challenging to derive. Our work commences by establishing two fundamental equivalence criteria pertaining to invariance properties. Building upon these, we devise two single-step algorithms grounded in LMI frameworks for the computation of RPI sets, denominated as Algorithm 1 and Algorithm 2. A comparative assessment between these algorithms is subsequently undertaken.

Experimental evaluations disclose that although both algorithms converge upon identical RPI sets, Algorithm 1 distinguishes itself through reduced optimization variable requirements and exhibits superior computational efficiency relative to Algorithm 2. This distinction carries substantial weight in practical implementations, as augmented computational speed directly contributes to enhanced real-time system performance and cost efficiency.

Capitalizing on the merits of Algorithm 1, we further propose an iterative refinement strategy aimed at incrementally approximating the target volume of RPI sets over successive iterations. This iterative algorithm ingeniously retains the computational efficiency hallmark of Algorithm 1 while simultaneously alleviating resultant conservatism, thereby striking a balance between precision and performance.

The structure of this paper is as follows: Preliminaries and problem formulation are given in Section $$ {\text{ 2.1 }} $$. Section $$ {\text{ 2.2 }} $$ presents data-based invariances and constraints. Section $$ {\text{ 2.3 }} $$ gives the relevant algorithms to maximize the RPI $$ {\text{sets}} $$. Numerical examples are given in Section $$ {\text{ 3 }} $$. Section 4 provides the discussion. Section 5 gives a summary of this article.

Notations: The following notations are made in this paper. The set of positive real numbers is denoted by $$ R_{+} $$; a diagonal matrix with positive members is denoted by $$ D_+^n $$; the matrix $$ A $$'s transpose is $$ A^T $$. $$ R^{n} $$ denotes $$ n $$-dimensional Euclidean space, and $$ R^{n \times m} $$ denotes a set of $$ n \times m $$-dimensional matrices. A matrix containing zeros in the relevant dimension is written as $$ \mathbf{0} $$. $$ \mathbf{1}_m $$ represents the vector of ones of dimension $$ m $$. The identity matrix in $$ m $$ dimensions is written as $$ I_m $$, and the $$ i $$-th column of the identity matrix is denoted by the symbol $$ e_i $$. $$ \ast $$'s represents the matrix element that is uniquely identifiable by symmetry. For the symmetry matrix $$ X $$, $$ X $$$$ \succeq $$ 0 ($$ \succ $$ 0) indicates that the matrix $$ X $$ is a semi-positive definite matrix (positive definite matrix). Let $$ A\in R^{m \times n} $$ be a matrix of $$ n $$ vectors $$ A $$ = [$$ a_1\dots a_n $$], we define the vectorization of $$ A $$ as $$ \vec{A} = [a_{1}^{T}\dots a_{n}^{T}]^T \in R^{mn} $$. For the finite set $$ \Theta_{v} $$ = $$ \left\lbrace \theta^1,\dots,\theta^t\right\rbrace $$ with $$ \theta ^j\in R^n $$ for $$ j=1,\dots, t $$, $$ conv(\Theta_{v}) = \left\lbrace \theta \in R^{n}:\theta =\sum_{j=1}^{t}\alpha_{j}\theta^{j},s.t \sum_{j=1}^{t}\alpha_{j}=1,\alpha_{j}\in [0,1]\right\rbrace $$ denotes the convex-hull of a $$ \Theta_{v} $$, and $$ A \otimes B $$ denotes the Kronecker product between $$ A $$ and $$ B $$.

2. METHODS

2.1. Preliminaries and problem formulation

2.1.1. Preliminaries

Lemma 1 (^[16]): (Vectorization) For matrices $$ A\in R^{m\times n} $$, $$ B\in R^{n\times l} $$, $$ C\in R^{l\times k} $$ and $$ D\in R^{m\times k} $$, the matrix equation $$ ABC = D $$ is equivalent to

(1a)

$$ \begin{equation} \left(C^{T}\otimes A\right) \vec{B} = \overrightarrow{ABC}, \end{equation} $$

(1b)

$$ \begin{equation} \overrightarrow{ABC} = \left( C^{T}B^{T}\otimes I_{k}\right)\vec{A}. \end{equation} $$

Lemma 2 (^[17]): (Schur complement): Given the matrix $$ S=\left[\begin{array}{cc} S_{1} & S_{2} \\ S_{3} & S_{4} \\ \end{array}\right] $$, where $$ S_{1} $$ is a positive definite matrix, and define the Schur complement matrix of $$ S_{1} $$ as $$ M = S_{4}-S_{3}S_{1}^{-1}S_{2} $$, then $$ S \succeq 0 (\succ 0) \Leftrightarrow M \succeq 0 (\succ 0) $$.

Definition 1 (^[1]): (Polyhedral set) A convex polyhedral set is a set of the form $$ \mathcal{P}(F,g) =\left\lbrace x :Fx\leq g \right\rbrace $$. A polyhedral set includes the origin if and only if $$ g \geq 0 $$ and includes the origin as an interior point if and only if $$ g>0 $$.

Definition 2 (^[1]): (0–Symmetric convex polyhedral set) A 0–symmetric convex polyhedral set can always be represented in the form $$ \mathcal{\bar{P}}(F,g) =\left\lbrace x :-g \leq Fx\leq g\right\rbrace $$. If $$ \mathcal{\bar{P}}(F,g) $$ includes $$ 0 $$ as an interior point, up to a normalization, it can be represented as $$ \mathcal{\bar{P}}(F,\mathbf{\bar{1}}) =\left\lbrace x :-\mathbf{\bar{1}} \leq Fx\leq \mathbf{\bar{1}} \right\rbrace $$, where $$ \mathbf{\bar{1}}^{T} =\left[1\enspace1 \dots 1 \right] $$.

$$ {\text{The properties of 0-symmetric convex polyhedral are as follows:}} $$

(ⅰ) If $$ { x \in \mathcal{P} , \text{then} -x \in \mathcal{P} .} $$

(ⅱ) $$ {\text{If} x,y \in \mathcal{P} , \text{then for any} \lambda \in [0,1], \lambda x+(1-\lambda)y \in \mathcal{P} } $$.

2.1.2. Problem formulation

The following linear discrete-time system $$ {\text{with no control inputs}} $$

(2)

$$ \begin{equation} x(k+1)=Ax(k)+w(k) \end{equation} $$

is considered in this paper, where $$ x(k)\in R^{n} $$$$ {\text{and}} $$$$ w(k) \in R^{n} $$ are the system state and the external disturbance at time $$ k $$, respectively. The $$ {\text{matrix}} $$$$ A $$ is unknown. A state trajectory of $$ T+1 $$ samples $$ \left\lbrace x(k)\right\rbrace _{k=1}^{T+1} $$ is generated from the system $$ \left(2 \right) $$. The generated data is denoted by:

(3a)

$$ \begin{aligned} X^{+} &\triangleq \left[ x(2)\enspace x(3)\dots x(T+1)\right] \in R^{n\times T}, \end{aligned} $$

(3b)

$$ \begin{aligned} X &\triangleq \left[x(1)\enspace x(2) \dots x(T)\right] \in R^{n\times T}. \end{aligned} $$

The constraints of the system state and disturbances are as follows

(4a)

$$ \begin{aligned} \mathcal{X}& \triangleq\left\lbrace x\in R^{n}:Fx\leq \mathbf{1}_{n_x}\right\rbrace, \end{aligned} $$

(4b)

$$ \begin{aligned} \mathcal{W}& \triangleq\left\lbrace w\in R^{n}: \mid Dw \mid\leq \mathbf{1}_{n_w}\right\rbrace, \end{aligned} $$

where $$ F \in R^{n_{x} \times n} $$ and $$ D \in R^{n_{w} \times n} $$ are known. Note that the external disturbances of the system are bounded.

Define a set of feasible models as follows:

(5)

$$ \begin{equation} \begin{split} \mathcal{A}\triangleq\left\lbrace A\in R^{n\times n}:x(k+1)-Ax(k) \in \mathcal{W}, k=1,\dots,T\right\rbrace \end{split} \end{equation} $$

By using the $$ {\text{matrices}} $$ defined in (3) and the disturbance set $$ \mathcal{W} $$ in (4), (5) can be expressed as follows:

(6)

$$ \begin{equation} \mathcal{A}\triangleq\left\lbrace A\in R^{n\times n}:-\mathbf{\bar{1}}\leq DX^{+}-DAX\leq \mathbf{\bar{1}} \in \mathcal{W}\right\rbrace, \end{equation} $$

where $$ \mathbf{\bar{1}} \triangleq\left[\mathbf{1}_{n_w} \enspace \mathbf{1}_{n_w} \enspace\dots\mathbf{1}_{n_w} \right] \in R^{n_{w}\times T} $$. $$ {\text{By using}} $$ Lemma 1, (6) is rewritten as follows,

(7)

$$ \begin{equation} \mathcal{A}\triangleq\left\lbrace A\in R^{n\times n}:-\mathbf{\bar{1}}_{Tn_{w}}+d\leq Z\vec{A}\leq \mathbf{\bar{1}}_{Tn_{w}}+d \in \mathcal{W}\right\rbrace, \end{equation} $$

where

(8a)

$$ \begin{aligned} Z&\triangleq\left(X^{T}\otimes D\right) \in R^{T_{n_w}\times n^{2}}, \end{aligned} $$

(8b)

$$ \begin{aligned} d&\triangleq \left[ \begin{array}{c} Dx(2)\\ Dx(3)\\ \vdots\\ Dx(T+1) \end{array}\right] \in R^{T_{n_{w}}}. \end{aligned} $$

Proposition 1 ^{[12, 18]}: The set $$ \mathcal{A} $$ in $$ \left(6 \right) $$ is a bounded polyhedron if and only if $$ rank\left(X\right) =n $$ and $$ D $$ has a full column rank.

Remark 1: If the condition is not met, then set $$ \mathcal{A} $$ is unbounded, making it difficult to find a feasible RPI set.

$$ {\text{To enhance clarity, denote} \; x(k+1) \; \text{as} \; x^{+} , \text{then system (2) can be represented by}:} $$

(9)

$$ \begin{equation} x^{+}=Ax+w. \end{equation} $$

Consider the following set of 0–symmetric convex polyhedral set with predefined complexity $$ n_{p} $$, which is given as follows:

(10)

$$ \begin{equation} \mathcal{P}\triangleq\left\lbrace x \in R^{n}:-\mathbf{\bar{1}}\leq PW^{-1}x \leq \mathbf{\bar{1}}\right\rbrace, \end{equation} $$

where $$ P\in R^{n_{p}\times n} $$, $$ W\in R^{n \times n} $$. Assuming that $$ W $$ is reversible, this will be ensured by the invariance conditions in the form of LMIs.

The set $$ \mathcal{P} $$ is the RPI set of system (9), if the following condition is satisfied:

(11)

$$ \begin{equation} x\in \mathcal{P} \Rightarrow x^{+} \in \mathcal{P},\forall w \in\mathcal{W},\forall A \in \mathcal{A}. \end{equation} $$

$$ {\text{The set} \; \mathcal{P} \; \text{also must adhere to the state constraints, meaning that} \; \mathcal{P} \subseteq \mathcal{X} , \text{which leads to the following}:} $$

(12)

$$ \begin{equation} x\in \mathcal{P} \Rightarrow x\in\mathcal{X}. \end{equation} $$

$$ {\text{From the above (3)-(4), (9)-(12), the problem addressed in this article is formulated as follows:}} $$

Problem: Given the data in (3), state constraints in (4) and matrix $$ P $$, compute $$ W $$ in (10) such that:

(1) The invariance condition in (11) is satisfied;

(2) The set $$ \mathcal{P} $$ satisfies state constraints in (4);

(3) Maximize the volume of the set $$ \mathcal{P} $$.

2.2. Invariance conditions and constraints based on data

To render the state constraints of the system and invariance conditions more manageable, consider the following coordinate transformation^[14]:

(13)

$$ \begin{equation} \theta=W^{-1}x \Leftrightarrow x=W\theta. \end{equation} $$

With the coordinate transformation, (10) can be expressed as:

(14)

$$ \begin{equation} \mathcal{P}\triangleq\left\lbrace W\theta \in R^{n}:\theta \in \Theta\right\rbrace , \end{equation} $$

$$ {\text{where Theta is a set of 0-symmetric convex polyhedra defined as follows:}} $$

(15)

$$ \begin{equation} \Theta \triangleq\left\lbrace \theta \in R^{n}:-\mathbf{1}_{n_p}\leq P\theta \leq \mathbf{1}_{n_p} \right\rbrace. \end{equation} $$

For a fixed $$ P $$, the vertices of $$ \Theta $$ are known, so $$ \Theta $$ can be written as a convex hull of these finite vertices ^[1].

(16)

$$ \begin{equation} \Theta=conv\left( \left\lbrace \theta^{1},\dots,\theta^{2\sigma}\right\rbrace \right), \end{equation} $$

where $$ \theta^{1},\dots,\theta^{2\sigma} $$ are the vertices of the $$ \Theta $$, and $$ \sigma $$ is known and determined by matrix $$ P $$.

Due to the symmetry of the set $$ \Theta $$, there is $$ \theta^{j+\sigma} = -\theta^{j} $$, for $$ j = 1,\dots,\sigma $$.} And by using (13), the state constraints are rewritten in $$ \theta $$ state-space as follows

(17)

$$ \begin{equation} \begin{split} HW\theta\leq\mathbf{1}_{n_x}, \forall \theta \in \Theta\Leftrightarrow \\-\mathbf{1}_{n_x} \leq HW\theta^{j}\leq\mathbf{1}_{n_x}, j=1,\dots,\sigma. \end{split} \end{equation} $$

With (13), system (9) can be expressed as:

(18)

$$ \begin{equation} W\theta^{+}=AW\theta+w, \end{equation} $$

where $$ A\in \mathcal{A} $$, $$ w\in \mathcal{W} $$.

Two equivalent invariance conditions for system (18) in the $$ \theta $$ state space are shown below.

Lemma 3 (^[14]): For system $$ \left(18 \right) $$, where the set $$ \Theta $$ in (15) is RPI set, then the following two conditions are equivalent:

(ⅰ) for all $$ \theta \in \Theta $$, $$ \forall w\in \mathcal{W} $$ and $$ \forall A\in \mathcal{A} $$,

(19)

$$ \begin{equation} \theta^{+}=\left( W^{-1}AW\theta+W^{-1}w\right) \in \Theta \end{equation} $$

(ⅱ) for each vertex $$ \theta^{j} $$, $$ j=1,\dots,2\sigma $$ of the set $$ \Theta $$, $$ \forall w\in \mathcal{W} $$ and $$ \forall A\in \mathcal{A} $$,

(20)

$$ \begin{equation} \left( \theta^{j}\right) ^{+}=\left( W^{-1}AW\theta^{j}+W^{-1}w\right) \in \Theta \end{equation} $$

2.3. Maximize the volume of the RPI set

In this section, we propose data-driven sufficient LMI conditions for computing the matrix $$ W $$, ensuring that the set $$ \Theta $$ remains invariant.

First, we will consider condition (19) for RPI of the set $$ \Theta $$. Applying (1b), system (18) can be further written as follows:

(21)

$$ \begin{equation} W\theta^{+}=\underbrace{((W\theta)^{T}\otimes I_{n})}_{g(W,\theta)}\vec{A}+w. \end{equation} $$

To achieve fewer conservative LMI conditions, the variables $$ V_{i} \in R^{n\times n} $$ are introduced and signals $$ \xi_{i}=V_{i}^{-1}W(\theta)^{+} $$, for $$ i=1,\dots,n_{p} $$. The dynamics (21) can then be expressed as follows

(22)

$$ \begin{equation} g(W,\theta)\vec{A}+w-V_{i}\xi_{i}=\textbf{0}. \end{equation} $$

From (15), the invariance condition (19) is given as follows: for all $$ \theta \in \Theta $$,

(23)

$$ \begin{equation} 1-(e_{i}^{T}P\theta^{+})^{2} \geq 0,\forall w\in \mathcal{W}, \forall A\in \mathcal{A}, \end{equation} $$

where $$ e_{i} $$ is the identity matrix $$ I_{n_{p}} $$'s $$ i $$-th column vector. By using $$ \xi_{i}=V_{i}^{-1}W(\theta)^{+} $$, (23) can be written as follows

(24)

$$ \begin{equation} 1-(e_{i}^{T}PW^{-1}V_{i}\xi_{i})^{2} \geq 0,\forall w\in \mathcal{W}, \forall A\in\mathcal{A}. \end{equation} $$

We multiply (24) by positive scalar variable $$ \phi_{i}>0 $$ and lower bound it by terms that are known to be non-negative (S-procedure ^[15]). This gives

(25)

$$ \begin{equation} \begin{split} \phi_{i}(1-(e_{i}^{T}PW^{-1}V_{i}\xi_{i})^{2})\geq \underbrace{2\xi_{i}^{T}(g(W,\theta)\vec{A}+w-V_{i}\xi_{i})}_{0}\\ +\underbrace{((\textbf{1}+d)-Z\vec{A})^{T}\Lambda_{i}((\textbf{1}+d)+Z\vec{A})}_{\geq 0}\\+\underbrace{(\textbf{1}+Dw)^{T}\Gamma_{i}(\textbf{1}-Dw)}_{\geq 0} \end{split} \end{equation} $$

where $$ \Lambda_{i}\in D_{+}^{Tn_{w}} $$, $$ \Gamma_{i} \in D_{+}^{n_{w}} $$. In this way, we get the sufficient condition of (23), i.e., (25).

Through (22), (4), and (7), it can be proved that the right side of (25) is non-negative for all $$ w\in \mathcal{W} $$, for all $$ A\in \mathcal{A} $$. After that, (25) can be rewritten as the following quadratic form:

(26)

$$ \begin{equation} \textbf{y}^{T}T_{i}(W,{\phi_{i}},\Gamma_{i},\Lambda_{i},V_{i})\textbf{y} \succeq 0, \forall \textbf{y}, \end{equation} $$

where $$ \textbf{y}=[1\; \vec{A}^{T}\; w^{T}\; -\xi_{i}^{T}]^{T} $$ and matrix $$ T_{i} $$ is symmetrical. The invariance condition holds, if $$ T_{i} $$$$ \succeq 0 $$; i.e.,

(27)

$$ \begin{aligned} T_{i}\triangleq\left[ \begin{array}{ccccc} r_{i}&-d^{T}\Lambda_{i}Z& \mathbf{0}&\mathbf{0}\\ *& Z^{T}\Lambda_{i }Z &0& g^{T}\left(W,\theta\right) \\ *&*&D^{T}\Gamma_{i}D&I_{n}\\ *&*&*&V_{i}+V_{i}^{T}-V_{i}^{T}L_{i}V_{i}\\ \end{array}\right] \succeq 0, \end{aligned} $$

where $$ L_{i}\triangleq\phi_{i}W^{-T}P^{T}e_{i}e_{i}^{T}PW^{-1} $$, $$ r_{i}\triangleq\phi_{i}-\textbf{1}^{T}\Lambda_{i}\textbf{1}-\textbf{1}_{n_{w}}^{T}\Gamma_{i}\textbf{1}_{n_{w}}+d^{T}\Lambda_{i}d $$ and $$ g(W,\theta) $$ is as defined in (21). Take note of the nonlinearity of the block(2, 4) in (27) and the nonlinear dependence of the block(4, 4) on $$ \phi_{i} ,V_{i} $$ and $$ W $$ in (27), which will be resolved by introducing new matrix variables and appropriate scaling.

Theorem 1: Given the matrices $$ X $$, $$ X^{+} $$ and $$ P $$, if there exists $$ W\in R^{n\times n} $$, $$ W_{1}\in R^{n} $$ and variables $$ \phi_{i} \in R_{+} $$, $$ \Lambda_{i}\in D_{+}^{Tn_{w}} $$, $$ \Gamma_{i} \in D_{+}^{n_{w}} $$, $$ X_{i}, V_{i} \in R^{n\times n} $$ that satisfy LMIs for $$ i=1,\dots n_{p} $$,

(28)

$$ \begin{equation} \left[ \begin{array}{cc} W+W^{T}-X_{i}&r_{i} P^{T}e_{i} \\ r_{i} e_{i}^{T} P&r_{i} \end{array}\right] \succeq 0, \end{equation} $$

(29)

$$ \begin{equation} \left[ \begin{array}{ccccc} r_{i}&-d^{T}\Lambda_{i}Z& \mathbf{0}&\mathbf{0}&\mathbf{0}\\ *& Z^{T}\Lambda_{i}Z &0& g^{T}\left(W_{1}\right)&\mathbf{0} \\ *&*&D^{T}\Gamma_{i}D&I_{n}&\mathbf{0}\\ *&*&*&V_{i}+V_{i}^{T}&V_{i}^{T}\\ *&*&*&*&X_{i} \end{array}\right] \succeq 0, \end{equation} $$

where

(30a)

$$ \begin{aligned} g(W_{1})&\triangleq\left( W_{1}\right) ^{T}\otimes I_{n} \in R^{n\times n^{2}}, \end{aligned} $$

(30b)

$$ \begin{aligned} W_{1}&\triangleq W\theta, \end{aligned} $$

(30c)

$$ \begin{aligned} r_{i}&\triangleq\phi_{i}-\mathbf{1}^{T}\Lambda_{i}\mathbf{1}-\mathbf{1}_{n_{w}}^{T}\Gamma_{i}\mathbf{1}_{n_{w}}+d^{T}\Gamma_{i}d \in R, \end{aligned} $$

then the set $$ \mathcal{P} $$ in $$ (14) $$ is RPI.

proof:

In order to resolve the nonlinearity in the block (2, 4) of (27), let us introduce a new matrix variable $$ W_{1}=W\theta $$ such that $$ g(W_{1}) $$ is linear.

First, let us prove the LMI condition (28) stated in Theorem 1. $$ {\text{Introduce}} $$ the positive-definite matrix variable $$ X_{i} $$ that satisfies

(31)

$$ \begin{equation} X_{i}^{-1}-L_{i} \succ 0 \Leftrightarrow X_{i}^{-1}-\phi_{i}W^{-T}P^{T}e_{i}e_{i}^{T}PW^{-1} \succ 0. \end{equation} $$

Applying Lemma 2 to (31), we obtain

(32)

$$ \begin{equation} \left[ \begin{array}{cc} X_{i}^{-1}&\phi_{i}W^{-T}P^{T}e_{i}\\ \phi_{i}e_{i}^{T}PW^{-1}&\phi_{i} \end{array}\right] \succ 0, \end{equation} $$

By applying congruence transformation to (32) using the congruence transformation matrix $$ diag\left\lbrace W,I_{n}\right\rbrace $$ which is full rank real matrix, then we get

(33)

$$ \begin{equation} \left[ \begin{array}{cc} W^{T}X_{i}^{-1}W&\phi_{i}P^{T}e_{i}\\ \phi_{i}e_{i}^{T}P&\phi_{i} \end{array}\right] \succ 0. \end{equation} $$

To address the nonlinearity in the block (1, 1) of (33), we perform the following steps:

(34)

$$ \begin{equation} \begin{split} W^{T}X_{i}^{-1}W &= {(W-X_{i})^{T}X_{i}^{-1}(W-X_{i})}\\&+W+W^{T}-X_{i}\\ &\succeq W+W^{T}-X_{i}, \end{split} \end{equation} $$

then the $$ W^{T}X_{i}^{-1}W $$ in (33) can be replaced by $$ W+W^{T}-X_{i} $$, and we obtain a sufficient LMI condition for (33) as given in (28).

Next, let us prove the LMI condition (29) stated in Theorem 1. From (31), the condition (27) can be rewritten as

35

$$ \left[\begin{array}{cccc} r_i & -d^T \Lambda_i Z & \mathbf{0} & \mathbf{0} \\ * & Z^T \Lambda_i Z & 0 & g^T\left(W_1\right) \\ * & * & D^T \Gamma_i D & I_n \\ * & * & * & V_i+V_i^T-V_i^T X_i^{-1} V_i \end{array}\right] \geq 0, $$

then using the Schur complement, we obtain (29).

Considering the condition (20) of the robust invariance of the set $$ \Theta $$, the following Theorem 2 can be obtained.

Theorem 2 (^[14]): Given the matrices $$ X $$, $$ X^{+} $$ and $$ P $$, if there exists $$ W\in R^{n\times n} $$ and variables $$ \phi_{ij}\in R_{+} $$, $$ \Lambda_{ij}\in D_{+}^{Tn_{w}} $$, $$ \Gamma_{ij} \in D_{+}^{n_{w}} $$, $$ X_{ij}, V_{ij} \in R^{n\times n} $$ that satisfy LMIs for $$ i=1,\dots n_{p} $$ and $$ j=1,\dots 2\sigma $$,

(36)

$$ \begin{equation} \left[ \begin{array}{cc} W+W^{T}-X_{ij}&r_{ij} P^{T}e_{i} \\ r_{ij} e_{i}^{T} P&r_{ij} \end{array}\right] \succeq 0, \end{equation} $$

(37)

$$ \begin{equation} \left[ \begin{array}{ccccc} r_{ij}&-d^{T}\Lambda_{ij}Z& \mathbf{0}&\mathbf{0}&\mathbf{0}\\ *& Z^{T}\Lambda_{ij }Z &0& g^{T}\left(W,\theta^{j}\right)&\mathbf{0} \\ *&*&D^{T}\Gamma_{ij}D&I_{n}&\mathbf{0}\\ *&*&*&V_{ij}+V_{ij}^{T}&V_{ij}^{T}\\ *&*&*&*&X_{ij} \end{array}\right] \succeq 0, \end{equation} $$

where,

(38a)

$$ \begin{aligned} g(W,\theta^{j})&\triangleq\left( W\theta^{j}\right) ^{T}\otimes I_{n} \in R^{n\times n^{2}}, \end{aligned} $$

(38b)

$$ \begin{aligned} r_{ij}&\triangleq\phi_{ij}-\mathbf{1}^{T}\Lambda_{ij}\mathbf{1}-\mathbf{1}_{n_{w}}^{T}\Gamma_{ij}\mathbf{1}_{n_{w}}+d^{T}\Gamma_{ij}d \in R, \end{aligned} $$

then the set $$ \mathcal{P} $$ in $$ (14) $$ is RPI.

We aim to identify the largest set $$ \mathcal{P} $$in (14) that satisfies the state constraints outlined in (4a) and the sufficient conditions for invariance specified in (28) and (29) (or alternatively, (36) and (37)). Given a matrix $$ P $$, it is known that the volume of $$ P $$ is proportional to its determinant, denoted as $$ \left|det(W) \right| $$^[19]. Consequently, we can determine the largest set $$ P $$ by formulating a SDP problem. Algorithms 1 and 2 based on the one-step method are given below.

Algorithm 1: Computing

$$ RPI $$

set.

Input:

$$ F $$

,

$$ D $$

,

$$ P $$

;
Input: Optimal values for

$$ \mathbf{W,X_{i}} $$

;
Objective function

$$ \underset{Z_{SDP}}{max}\enspace\enspace log\; det(W) $$

Optimization variables

$$ Z_{SDP}\triangleq \left(W,W_{1},X_{i},V_{i},\phi_{i},\Lambda_{i},\Gamma_{i}\right) $$

Symmetry constraint

$$ W=W' $$

(in order for the objective function to be concave)
State constraints (17)
Invariance conditions (28), (29)

Algorithm 2: Computing

$$ RPI $$

set.

Input:

$$ F $$

,

$$ D $$

,

$$ P $$

;
Output: Optimal values for

$$ \mathbf{W,X_{ij}} $$

;
Objective function

$$ \underset{Z_{SDP}}{max}\enspace\enspace log\; det(W) $$

Optimization variables

$$ Z_{SDP}\triangleq \left(W,X_{ij},V_{ij},\phi_{ij},\Lambda_{ij},\Gamma_{ij}\right) $$

Symmetry constraint

$$ W=W' $$

(in order for the objective function to be concave)
State constraints (17)
Invariance conditions (36), (37)

Therefore, in order to obtain the desired large volume of RPI $$ {\text{sets}} $$, we need to solve the determinant maximization problem under LMI conditions. However, this problem is easy to solve only if $$ W $$ is symmetric ^[20]. The symmetry of $$ W $$ would introduce conservatism ^[21]; thus, we introduce an $$ {\text{iterative algorithm}} $$ based on Algorithm 1. In an $$ {\text{iterative algorithm}} $$, $$ W $$ does not need to maintain symmetry, and the algorithm also reduces the conservatism caused by the introduction of (34) ^[19]. At the $$ k $$-th iteration, let $$ W_{k} $$ and $$ X_{ik} $$ be the values of the variables $$ W $$, $$ X_{i} $$. At each subsequent iteration, the volume of the RPI $$ {\text{set increases}} $$, i.e., $$ |det(W_{k+1})| \geq |det(W_{k})| $$, if the following LMI condition is imposed,

(39)

$$ \begin{equation} W^{T}W_{k}+W_{k}^{T}W-W_{k}^{T}W_{k} \succeq Z \succ 0, \end{equation} $$

where $$ Z=Z^{T} \in R^{n \times n} $$ is the new matrix variable.

From ^[22], we get

(40)

$$ \begin{equation} \begin{split} (W-X_{i}Z_{ik})^{T}X_{i}^{-1}(W-X_{i}Z_{ik})=W^{T}X_{i}^{-1}W\\-W^{T}Z_{ik}-Z_{ik}^{T}W+Z_{ik}^{T}X_{i}Z_{ik} \succeq 0, \end{split} \end{equation} $$

thus we can obtain:

(41)

$$ \begin{equation} W^{T}X_{i}^{-1}W \succeq W^{T}Z_{ik}+Z_{ik}^{T}W-Z_{ik}^{T}X_{i}Z_{ik}, \end{equation} $$

where $$ Z_{ik} \triangleq X_{ik}^{-1}W_{k} $$, and then the linear term to the right of (41) can be used to substitute for the nonlinear term $$ W^{T}X_{i}^{-1}W $$ in (29). From (41), the condition (28) can be rewritten as follows,

(42)

$$ \begin{equation} \left[ \begin{array}{cc} W^{T}Z_{ik}+Z_{ik}^{T}W-Z_{ik}^{T}X_{i}Z_{ik}&r_{i}P^{T}e_{i} \\ {\text{ r_{i}e_{i}^{T}P }} &\phi_{i} \end{array}\right] \succeq 0 . \end{equation} $$

Therefore, we can obtain the following $$ {\text{iterative algorithm}} $$ based on Algorithm 1:

iterative algorithm Computing PRI set.
Input:$$ F $$, $$ D $$, $$ P $$, $$ W^{q} $$, $$ X_{i}^{q} $$; Output: Optimal values for $$ \mathbf{W} $$; Objective function $$ \underset{Z_{SDP}}{max}\enspace\enspace log\; det(Z) $$ Optimization variables $$ Z_{SDP}\triangleq \left(W,X_{i},V_{i},\phi_{i},\Lambda_{i},\Gamma_{i},Z \right) $$ State constraints (17) Invariance conditions (28), (42)

Remark 2: Algorithm 2 has $$ 10\sigma n_{p}+1 $$ optimization variables, the invariance conditions (36), (37) consist of $$ 2\sigma n_{p} $$ linear matrix inequalities respectively. The iterative algorithm has $$ 5n_{p}+2 $$ optimization variables; the invariance conditions (28), (42) consist of $$ n_{p} $$ linear matrix inequalities respectively. This suggests that both the iterative algorithm has fewer optimization variables and fewer LMIs for invariant conditions than Algorithm 2. In addition, Algorithm 1 is more conservative than iterative algorithm (see Example 2).

3. RESULTS

The algorithms in the experiments are implemented in Matlab by using CVX ^[23] and solver SeDuMi. And the MPT toolbox is used to manage polytopes ^[24]. The following liner discrete-time system is considered:

(43)

$$ \begin{equation} x(k+1)=\underbrace{\left[ \begin{array}{cc} 1&1 \\ 0&1 \end{array}\right]}_{A} x(k)+w(k). \end{equation} $$

The matrice $$ A $$ is unknown, but it is only used to gather data. $$ T $$ = 20 samples of data are collected by system (43), as shown in Figure 1. Assume the disturbance $$ w $$ in system (43) ranges between -0.1 and 0.1, with the state constraints being $$ (x_{1}, x_{2}) \in [-2, 2] \times [-2,2] $$.

Computation of robust positively invariant set based on direct data-driven approach

Figure 1. The data samples collected from the liner discrete-time system (43).

Example 1: By taking different matrices $$ P $$, the complexity of the RPI set is denoted as $$ n_{p} $$ = 2, 3, 4 respectively, and the corresponding matrices are computed as follows

$$ P_{2}=\left[ \begin{array}{cc} 1&0 \\ 0&1 \end{array}\right] $$, $$ P_{3}=\left[ \begin{array}{ccc} 10&10 \\ 10&0\\ 1&11 \end{array}\right] $$, $$ P_{4}=\left[ \begin{array}{cccc} -18&-55 \\ 18&55\\ 55&-18\\ 55&18 \end{array}\right] $$.

The matrice $$ W $$ is computed by running Algorithms 1 and 2 as follows:

$$ W_{21}=\left[ \begin{array}{cc} 3.5&0 \\ 0&1 \end{array}\right] $$(Algorithm 1, 2.039 $$ {\text{s}} $$), $$ W_{22}=\left[ \begin{array}{cc} 3.5&0 \\ 0&1 \end{array}\right] $$(Algorithm 2, 2.770 $$ {\text{s}} $$),

$$ W_{31}=\left[ \begin{array}{cc} 35&1 \\ 1&11 \end{array}\right] $$(Algorithm 1, 2.418 $$ {\text{s}} $$), $$ W_{32}=\left[ \begin{array}{cc} 35&1 \\ 1&11 \end{array}\right] $$(Algorithm 2, 4.552 $$ {\text{s}} $$),

$$ W_{41}=\left[ \begin{array}{cc} 151.4106&18.1135 \\ 18.1135&55.0114 \end{array}\right] $$(Algorithm 1, 2.410 $$ {\text{s}} $$), $$ W_{42}=\left[ \begin{array}{cc} 151.4106&18.1135 \\ 18.1135&55.0114 \end{array}\right] $$(Algorithm 2, 4.725 $$ {\text{s}} $$),

where $$ W_{21} $$ and $$ W_{22} $$ are obtained from $$ P_{2} $$; $$ W_{31} $$ and $$ W_{32} $$ are obtained by taking $$ P_{3} $$; $$ W_{41} $$ and $$ W_{42} $$ are obtained by taking $$ P_{4} $$.

It is discovered from experiment results that while Algorithm 2 yields the same results, Algorithm 1 runs faster and requires fewer optimization variables (see Remark $$ {\text{2}} $$) than Algorithm 2.

Example 2: The RPI sets with complexities $$ n_{p} $$ = 2, 3, 4 (i.e., $$ P_{2} $$, $$ P_{3} $$, $$ P_{4} $$ are selected respectively) which are obtained by Algorithm 1 are as follows

$$ W_{2}=\left[ \begin{array}{cc} 3.5&0 \\ 0&1 \end{array}\right] $$, $$ W_{3}=\left[ \begin{array}{cc} 35&1\\ 1&11 \end{array}\right] $$, $$ W_{4}=\left[ \begin{array}{cc} 151.4106&18.1135 \\ 18.1135& 55.0114 \end{array}\right] $$.

The corresponding RPI sets obtained by Algorithm 1 are shown in Figure 2.

Figure 2. Maximum volume RPI sets $$ \mathcal{P} $$ with different complexities: $$ n_{p} $$ = 2 (left), $$ n_{p} $$ = 3 (right), $$ n_{p} $$ = 4 (bottom).

After five times iterations of the $$ {\text{iterative algorithm}} $$, the RPI sets $$ \mathcal{P} $$ with complexities $$ n_{p} $$ = 2, 3, 4 are obtained as

$$ W_{2}=\left[ \begin{array}{cc} 3.5&0 \\ 0&1 \end{array}\right] $$, $$ W_{3}=\left[ \begin{array}{cc} 35&1\\ 1&11 \end{array}\right] $$, $$ W_{4}=\left[ \begin{array}{cc} 192.5402&-63.3969 \\ 18.1136& 55.0115 \end{array}\right] $$.

The corresponding RPI sets obtained by iterative algorithm are shown in Figure 3.

Figure 3. Maximum volume RPI sets $$ \mathcal{P} $$ with different complexities: $$ n_{p} $$ = 3 (left), $$ n_{p} $$ = 4 (right).

We find that the volume of RPI $$ {\text{sets}} $$ obtained by the $$ {\text{iterative algorithm }} $$ is larger than that obtained by Algorithm 1 (see Table 1 and Table 2), which indicates that Algorithm $$ 1 $$ is more conservative than the $$ {\text{iterative algorithm}} $$. For the given matrix $$ P $$, $$ {\text{the volume of} \mathcal{P} } $$ is proportional to the determinant $$ \left|det(W)\right| $$^[19]. And as $$ n_{p} $$ rises, the absolute value of the determinant of $$ W $$ grows; thus, the volume of the RPI set increases (see Table 1 and Figure 4), so that an invariant set with a higher volume can be obtained by adding $$ n_{p} $$ as a parameter.

Table 1

$$ \mathit{|det(W)|} $$ (Algorithm 1) $$ \mathit{vs.} $$ complexity of the matrix $$ P $$

complexity	$$ \mathit{\boldsymbol{n_{p}=2}} $$	$$ \mathit{\boldsymbol{n_{p}=3}} $$	$$ \mathit{\boldsymbol{n_{p}=4}} $$
$$ \|det(W)\| $$	3.5	384.0	8001.2

Table 2

$$ \mathit{|det(W)|} $$ (iterative algorithm) $$ \mathit{vs.} $$ complexity of the matrix $$ P $$

complexity	$$ \mathit{\boldsymbol{n_{p}=2}} $$	$$ \mathit{\boldsymbol{n_{p}=3}} $$	$$ \mathit{\boldsymbol{n_{p}=4}} $$
$$ \|det(W)\| $$	3.5	385.0	11740.3

Figure 4. $$ |det(W)|\; {vs.} $$ iterations of iterative algorithm.

4. DISCUSSION

A direct data-driven method is proposed to calculate the robust positive invariant (RPI) sets for discrete-time systems. To further reduce conservatism, an iterative algorithm based on Algorithm 1 is introduced. This approach does not require prior knowledge of the model or system identification. Future research could extend this methodology to nonlinear systems. Additionally, while the RPI sets discussed here are symmetric, the exploration of asymmetric RPI sets is another potential area for investigation.

5. CONCLUSIONS

This paper proposes a direct data-driven method to calculate RPI sets for discrete-time systems by deriving a set of invariance conditions expressed as linear matrix inequalities (LMIs). Subsequently, we maximize the volume of the invariant set using a SDP problem. We have developed two one-step algorithms based on LMIs to compute the RPI sets; experimental results indicate that Algorithm 1 requires fewer optimization variables compared to Algorithm 2 and demonstrates superior computational efficiency. Additionally, we have introduced an iterative approach based on Algorithm 1 to further reduce conservatism. Numerical examples have verified the effectiveness of the proposed algorithm.

DECLARATIONS

Acknowledgments

We would like to express our great appreciation to the editors and reviewers.

Authors' contributions

Conceptualization, methodology, resources, supervision: Yang H

Software, validation, writing - original draft preparation: Du Q

Formal analysis, investigation, writing - review and editing, visualization: Yang H, Du Q

All authors have read and agreed to the published version of the manuscript.

Availability of data and materials

Not applicable.

Financial support and sponsorship

None.

Conflicts of interest

Both authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

REFERENCES

1. Blanchini F, Miani S. Set-theoretic methods in control, 3th ed. Springer, 2008. Available from: https://link.springer.com/book/10.1007/978-3-319-17933-9[Last accessed on 28 Nov 2024].

2. Castelan EB, Tarbouriech S. Positively invariant polyhedral sets for discrete-time singular systems with additive perturbations. In Proceedings of the 35th IEEE Conference on Decision and Control; 13 December 1996; Kobe, Japan.

3. Blanchini F. Set invariance in control. Automatica 1999;35:1747-67.

4. Kolmanovsky I, Gilbert EG. Theory and computation of disturbance invariant sets for discrete‐time linear systems. Math Probl Eng 1998;4:317-67.

5. Raković SV, Kerrigan EC, Mayne DQ, Kouramas KI. Optimized robust control invariance for linear discrete-time systems: theoretical foundations. Math Probl Eng 2007;43:831-41.

6. Wang P, Feng X, Li W. Design of robust positively invariant set for nonholonomic vehicle. In Proceedings of the 2016 Chinese Control and Decision Conference (CCDC); 28-30 May 2016; Yinchuan, China.

7. Esterhuizen W, Aschenbruck T, Streif S. On maximal robust positively invariant sets in constrained nonlinear systems. Automatica 2020;119:109044.

8. Rakovic SV, Kerrigan EC, Kouramas K, Mayne DQ. Invariant approximations of the minimal robust positively invariant set. IEEE Trans Automat Contr 2005;50:406-10.

9. Trodden P. A one-step approach to computing a polytopic robust positively invariant set. IEEE Trans Automat Contr 2016;61:4100-5.

10. Xue B, Zhan N. Robust invariant sets computation for discrete-time perturbed nonlinear systems. IEEE Trans Automat Contr 2022;67:1053-60.

11. Hou ZS, Wang Z. From model-based control to data-driven control: survey, classification and perspective. Inf Sci 2013;235:3-35.

12. Zhong B, Zamani M, Caccamo M. Synthesizing safety controllers for uncertain linear systems: a direct data-driven approach. In Proceedings of the 35th IEEE Conference on Decision and Control; 23-25 August 2020; Trieste, Italy.

13. Bisoffi A, De Persis C, Tesi P. Data-based guarantees of set invariance properties. IFAC-PapersOnLine 2020;53:3953-8.

14. Mejari M, Gupta A. Direct data-driven computation of polytopic robust control invariant sets and state-feedback controllers. In Proceedings of the 62nd IEEE Conf. on Decision and Control (CDC); 13-15 December 2023; Singapore, Singapore.

15. Scherer CW. A full block S-procedure with applications. In Proceedings of the 36th IEEE Conf. on Decision and Control; 12 December 1997; San Diego, CA, USA.

16. Abadir KM, Magnus JR. Matrix algebra Cambridge University Press; 2005.

17. Cottle RW. Manifestations of the Schur complement. Linear Algebra Appl 1974;8:189-211.

18. Willems JC, Markovsky I, Rapisarda P, De Moor BLM. A note on persistency of excitation. Syst Control Lett 2005;54:325-9.

19. Gupta A, Köroğlu H, Falcone P. Computation of robust control invariant sets with predefined complexity for uncertain systems. Intl J Robust Nonlinear 2021;31:1674-88.

20. Vandenberghe L, Boyd S, Wu SP. Determinant maximization with linear matrix inequality constraints. SIAM J Matrix Anal Appl 1998;19:499-533.

21. Gupta A, Köroğlu H, Falcone P. Restricted-complexity characterization of control-invariant domains with application to lateral vehicle dynamics control. In Proceedings of the 2017 IEEE 56th Annual Conference on Decision and Control (CDC); 12-15 December 2017; Melbourne, VIC, Australia.

22. Gupta A, Köroğlu H, Falcone P. Computation of low-complexity control-invariant sets for systems with uncertain parameter dependence. Automatica 2019;101:330-7.

23. CVX: matlab software for disciplined convex programming; 2020. Available from: http://cvxr.com/cvx[Last accessed on 28 Nov 2024].

24. Herceg M, Kvasnica M, Jones CN, Morari M. Multi-parametric toolbox 3.0. In Proceedings of the 2013 European Control Conference (ECC); 17-19 July 2013; Zurich, Switzerland.

Cite This Article

Research Article

Open Access

Computation of robust positively invariant set based on direct data-driven approach

Qian Du, Hongli Yang

How to Cite

Du, Q.; Yang, H. Computation of robust positively invariant set based on direct data-driven approach. Complex Eng. Syst. 2024, 4, 24. http://dx.doi.org/10.20517/ces.2024.76

Download Citation

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click on download.

Export Citation File:

RIS BibTeX EndNote

Type of Import

Direct Import Indirect Import

Tips on Downloading Citation

This feature enables you to download the bibliographic information (also called citation data, header data, or metadata) for the articles on our site.

Citation Manager File Format

Use the radio buttons to choose how to format the bibliographic data you're harvesting. Several citation manager formats are available, including EndNote and BibTex.

Type of Import

If you have citation management software installed on your computer your Web browser should be able to import metadata directly into your reference database.

Direct Import: When the Direct Import option is selected (the default state), a dialogue box will give you the option to Save or Open the downloaded citation data. Choosing Open will either launch your citation manager or give you a choice of applications with which to use the metadata. The Save option saves the file locally for later use.

Indirect Import: When the Indirect Import option is selected, the metadata is displayed and may be copied and pasted as needed.

About This Article

Copyright

© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, sharing, adaptation, distribution and reproduction in any medium or format, for any purpose, even commercially, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Data & Comments

Data

Views

72

Downloads

15

Citations

0

Comments

0

2

Comments

Comments must be written in English. Spam, offensive content, impersonation, and private information will not be permitted. If any comment is reported and identified as inappropriate content by OAE staff, the comment will be removed without notice. If you have any queries or need any help, please contact us at support@oaepublish.com.