Unified CACSD Toolbox for Hybrid Simulation and Robust Controller Synthesis with Applications in DC-to-DC Power Converter Control

Abstract

The current article presents the design, implementation, validation, and use of a Computer-Aided Control System Design (CACSD) toolbox for nonlinear and hybrid system uncertainty modeling, simulation, and control using $\mu$ synthesis. Remarkable features include generalization of classical system interconnection operations to nonlinear and hybrid systems, automatic computation of equilibrium points for nonlinear systems, and optimization of least conservative uncertainty bounds, with direct applicability for $\mu$ synthesis. A unified approach is presented for the step-down (buck), step-up (boost), and single-ended primary-inductor (SEPIC) converters to showcase the use and flexibility of the toolbox. Robust controllers were computed by minimization of the ${\mathcal{H}}_{\infty }$ norm of the augmented performance systems, encompassing a wide range of uncertainty types, and have been designed using the well-known mixed-sensitivity closed loop shaping $\mu$ synthesis method.

1. Introduction

Robust Control represents a massive point of interest when it comes to Control Theory, which has been heavily studied over the past decades. However, albeit Robust Control brings many benefits, it is still an open field in research which gathers increasingly more approaches over time. Basically, the goal of a robust controller is to accomplish a specified set of performances for bounded model uncertainties which can occur in practice due to various reasons. In other words, closed loop stability and performance are maintained even for model parameter variations and unmodeled dynamics alike.

Over the years, multiple and various approaches for designing robust controllers have been presented, some of them being implemented into dedicated toolboxes, such as MATLAB’s Robust Control Toolbox [1]. This toolbox gathers the most efficient ones based on ${\mathcal{H}}_2$, ${\mathcal{H}}_{\infty }$, and $\mu$ synthesis methods, and it is often considered a reference in research. However, while using these types of toolboxes leads to controllers which are optimal for their prescribed criterion, they are not necessarily best in terms of conventional performances. Additionally, of great use for defining and optimizing difficult robust control problems is the Global Optimization Toolbox from in [2], providing ready-to-use solvers using various state-of-the-art algorithms, such as Particle Swarm Optimization (PSO) and Genetic Algorithms (GA). An important work in this direction, for computing optimal weighting functions for the generalized plant model, is presented in [3].

Even though there is a large variety of CACSD toolboxes in the field, their number is still expanding due to the necessity of overcoming drawbacks that the already existing ones have. At this point, the purpose of new toolboxes is not only to determine robust controllers for a specific process class, but to use a unified approach that would make them work for more types of systems, even multiple interconnected systems, in various configurations. User experience is also more accentuated, which is why some of them incorporate graphical user interfaces (GUIs), for improved usability.

An example is Multivar, which is a MATLAB-based application used for multiple-input and multiple-output (MIMO) control design, presented in [4]. This toolbox supports two working modes. It allows the user to work both in function and GUI mode (which represents a configuration wizard for determining the controller). Multivar can be used for LTI systems with or without time delay and it allows creating a model; converting, approximating, and analyzing it; input–output pairing and decoupling; and controller design and evaluation. Besides this, the user is able to export the control design and compare it with other saved designs. Another GUI-based robust controller design tool, which was created in LabVIEW, is presented in [5], based on the ${\mathcal{H}}_{\infty }$ loop shaping method. However, the goal was to provide a simple, user-friendly interface to make it easier to use, especially for educational purposes. Therefore, as mentioned by the authors, it does not provide the same flexibility as other design tools on the market.

LCToolbox, as presented in [6], is another MATLAB software package which is used for robust controller design. One of the advantages of using this toolbox instead of classic MATLAB routines is the fact that it gathers all necessary steps for controller design in one place, while cutting the need of preprocessing steps such as separate construction of the plant, and postprocessing steps, such as closed loop simulation. LCToolbox can be used for both linear time-invariant (LTI) and linear parameter-varying (LPV) models, and it also incorporates system identification methods. The controller is obtained by using the ${\mathcal{H}}_{\infty }$ loop shaping method. Other ${\mathcal{H}}_{\infty }$-based CACSD toolboxes have been presented over the past years. One example is represented in [7], which is based on linearizing or convexifying the conventional non-convex constraints on the classical robustness margins of ${\mathcal{H}}_{\infty }$ constraints. The controller parameters are then computed by using an optimization solver. This toolbox was created for MATLAB, and some of its main features are represented by the large variety of control problems in which it could be used, such as multi-model systems; the toolbox is designed to work with the output data of MATLAB’s System Identification Toolbox [8]. The output of the toolbox is represented by a PID controller, which can be easily implemented. Another example of a ${\mathcal{H}}_{\infty }$-based CACSD toolbox is shown in [9], in which the main advantage is the reduced conservatism of almost all types of model uncertainties which are defined.

Controller order is an important factor when implementing it on real systems. Therefore, this might be an issue in some cases. However, methods that are determining a fixed structure controller are already presented, such as in [10], which is based on the ${\mathcal{H}}_2$ controller design method, but can be cumbersome to compute. In order to deal with the high order controller problem, other toolboxes include controller simplification steps to avoid the necessity of postprocessing, as presented in [11].

Currently acknowledged problems in this domain regard closed loop simulation, where performance validation is generally treated ad hoc, from one control problem to another. Another difficulty encountered is when the test cases were done only on the linearized system for which the controller is designed, without checking if the initially proposed performance values are additionally verified for the nonlinear plant model, and, also, uncertainty modeling is a very cumbersome operation. The purpose of the paper is to provide means for treating the previously stated problems in a unified manner, such as implementing automated testing, performance validation, and report generation.

In this current iteration of the toolbox, robust controllers were designed using the well-established routines from in [1]. The interface is scalable and the control logic and validation can be replaced with other user-defined methods, or the current robust control approach can be replaced with open source alternatives for the ${\mathcal{H}}_2$ and ${\mathcal{H}}_{\infty }$ optimization problems, such as presented in the thesis [12], with the possibility to refine the necessary solutions of the Algebraic Riccati Equations (AREs) using the algorithms from [13], while for the $\mu$ synthesis problem, the thesis [14] provides a flexible, open-source implementation using linear matrix inequalities (LMIs). A clear advantage over the ARE approach is that LMIs are capable of solving singular and close to singular problems. Alternatively, a mixed ${\mathcal{H}}_2$/${\mathcal{H}}_{\infty }$ approach for stabilization and optimization using fixed-order controllers can be found in [15]. As such, the current iteration of the proposed toolbox is MATLAB-dependent for certain key functionalities, especially with regards to numerical simulation, robust control, and optimization, although the exposed ideas and mathematical framework can be directly implemented in other software environments, such as Python, Scilab, or LabVIEW.

The remainder of the paper is structured in the following manner. Section 2 describes the software structure and features of the proposed toolbox; Section 3 describes a proposed end-to-end workflow exemplified using modeling, control, validation, and simulation of several DC-to-DC converter topologies; and Section 4 illustrates comparative discussions, proposed improvements, and completions for future work and conclusions.

2. Toolbox Structure and Functionalities

The proposed toolbox has been designed with the target of end-to-end design and implementation of closed loop control systems, starting from the definition of the uncertainty set of plants to be controlled, their required operating point, along with control performance specifications and controller synthesis, and ending with the controller validation for the initial desired plant set.

2.1. Toolbox Features

Proposed features and advantages over existing toolboxes available in the literature:

• specify finite-dimensional dynamical systems with the general framework from Equation (1) to be used with the MATLAB ode framework; ability to interconnect such systems in series, parallel, and linear-fractional transformations; this functionality is described in Section 2.2.1 and Section 2.2.2;
• specify hybrid dynamical systems in the framework from in [16] as in Equation (4), with the ability to interconnect such systems in series, parallel, and linear-fractional transformations, upper and lower; this feature is described in Section 2.2.3;
• automatically compute equilibrium points numerically, with the possibility to impose certain states, inputs, and/or outputs, while the remaining ones are deduced through numerical optimization; this feature is presented in Section 2.4;
• automatically compute the uncertainty model as requested alongside a nominal plant: additive, inverse additive, input and output multiplicative, etc. using a global optimization algorithm, such as particle swarm optimization, to be directly used as necessary for robust synthesis methods; removes the burden for the control engineer to manually do this process for each plant; this feature is presented in Section 2.5;
• flexible and scalable, all features are implemented through MATLAB code and does not need the use of Simulink, which can become cumbersome when treating families of plants and not a single, specific, plant at a time; also, to account for the operating point in the case of linearized, nonlinear, and hybrid systems, alike, the same interface for Model-in-the-Loop simulation is provided in the toolbox, as shown in Section 2.2 and Section 2.3;
• besides the automatic validation of the frequency response for the desired operating point of the linearized plant family, the toolbox runs tests accounting for the uncertainty behavior of the desired nonlinear plant, not only on the linearization which the controller has been designed for. Every specification imposed in the designed phase will be automatically tested for the entire nonlinear system family, as illustrated in the case studies from Section 3.

2.2. Systems Specification

The scope of software classes implemented and described in this section aims to provide a flexible framework for simulation by using the ordinary-differential equation ode solver exclusively, with the low-level requirement of integrating a differential state equation. As such, exogenous signals would be reference signals and disturbances, known a priori in a simulation context. The intrinsic signals, i.e., commands and corresponding measurements, are passed to their corresponding subsystems by means of ode. Figure 1 encompasses an overview of the toolbox classes described in Section 2.2, Section 2.3 and Section 2.4. When the relationship between two classes is of type inheritance, the inherited class will not redundantly recall all previous properties and methods from the base class in the diagram, unless they overload the methods and is explicitly noted.

2.2.1. Nonlinear Systems

For the purpose of this paper, we will focus on finite-dimensional systems: deterministic and stochastic. The so-called explicit or standard system form is obtained by writing the plant model in the following canonical form, using a set of differential equations and a set of output equations:

$$\{\begin{array}{}\mathrm{(1a)}& \dot{\mathbf{x}}\left(t\right)=F\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),t\right);\hfill \mathrm{(1b)}& \mathbf{y}\left(t\right)=h\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),t\right),\hfill \end{array}$$

with the vector maps F and h being Lipschitz functions. The input signal $\mathbf{u}\left(t\right)$ has dimension m, state signal $\mathbf{x}\left(t\right)$ with dimension n, and output signal $\mathbf{y}\left(t\right)$ with dimension p, with $t\ge 0$. The initial conditions of the system are $\mathbf{x}\left(0\right)={\mathbf{x}}_0\in {\mathbb{R}}^n$.

Dynamical systems of the form (1) are implemented in class System. This will be the baseline interface for all systems the toolbox works with. Its most important methods are sim, findEqPoint, and linearize, which will be briefly described. The method sim simulates the dynamical system described by Equation (1a) from the initial condition x0, using the exogenous signal u(t), which is a predetermined anonymous function with at least the input argument time. tfin can be a scalar time value representing the final simulation time, a simulation interval, or a vector of predetermined time values. The solver options and type are based on MATLAB’s ode framework options and are sent directly to it. The solver type can be selected from any of the supported functions: ode113, ode15s, ode15i, ode23, ode23t, ode23s, ode23tb, or ode45. After integrating the state equation, the output signal $\mathbf{y}\left(t\right)$ can be directly computed using the memoryless function h from (1b). A useful particularization is also the method simInitCond, with the only difference being that it replaces the time-varying input signal $\mathbf{u}\left(t\right)$ with a constant value ${\mathbf{u}}_0$, thus obtaining an impulse response. The method findEqPoint deduces an equilibrium point for the system given a set of specifications on the input, state, or output vectors and is described in detail in Section 2.4. After obtaining a valid equilibrium point, a linearized system can be obtained using the method linearize, also described there.

2.2.2. Linear Systems

Of particular interest for the framework and for control systems in general are linear time-invariant systems, which inherit the software interface from the System class, are implemented in the class LTISystem and are defined by

$$\{\begin{array}{}\mathrm{(2a)}& \dot{\mathbf{x}}\left(t\right)=A\mathbf{x}\left(t\right)+B\mathbf{u}\left(t\right);\hfill \mathrm{(2b)}& \mathbf{y}\left(t\right)=C\mathbf{x}\left(t\right)+D\mathbf{u}\left(t\right).\hfill \end{array}$$

Separately, a nonlinear system can be linearized in the vicinity of an operating point, which is an equilibrium point for said system. The operating point $\left({\mathbf{u}}_0,{\mathbf{x}}_0,{\mathbf{y}}_0,t_0\right)$ can be provided by the user or can be computed using the functionality from Section 2.4. The linearized system will work with variations of the initial variables and have the following model:  

$$\left\{\begin{array}{c}\Delta \dot{\mathbf{x}}\left(t\right)=A·\Delta \mathbf{x}\left(t\right)+B·\Delta \mathbf{u}\left(t\right);\hfill \\ \Delta \mathbf{y}\left(t\right)=C·\Delta \mathbf{x}\left(t\right)+D·\Delta \mathbf{u}\left(t\right);\hfill \end{array}\right.\iff \left\{\begin{array}{c}\dot{\mathbf{x}}\left(t\right)=A\left(\mathbf{x}\left(t\right)-{\mathbf{x}}_0\right)+B\left(\mathbf{u}\left(t\right)-{\mathbf{u}}_0\right);\hfill \\ \mathbf{y}\left(t\right)=C\left(\mathbf{x}\left(t\right)-{\mathbf{x}}_0\right)+D\left(\mathbf{u}\left(t\right)-{\mathbf{u}}_0\right)+{\mathbf{y}}_0.\hfill \end{array}\right.$$

(3)

This latter structure is useful for MiL simulations and is implemented in the auxiliary class LTIEqSystem, seen as an affine nonlinear system. The great advantage of having the system from Equation (3) readily available is that it is interchangeable with the initial nonlinear interface in a closed loop context without making further adaptations in the source code and can be used to study the performance degradation obtained by replacing the controller from the linearized system to the nonlinear plant.

2.2.3. Hybrid Systems

A useful extension of framework (1) for hybrid systems, to account for system discontinuities, is with structures described in [16,17]:

$$\{\begin{array}{}\mathrm{(4a)}& \dot{\mathbf{x}}\left(t\right)\hfill & =F\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),t\right),\left(\mathbf{x},\mathbf{u},t\right)\in \mathcal{C};\hfill \mathrm{(4b)}& {\mathbf{x}}^{+}\left(t\right)\hfill & =G\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),t\right),\left(\mathbf{x},\mathbf{u},t\right)\in \mathcal{D};\hfill \mathrm{(4c)}& \mathbf{y}\left(t\right)\hfill & =h\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),t\right),\hfill \end{array}$$

with $F:{\mathbb{R}}^{n+m+1}\to {\mathbb{R}}^n$ as the flow function, $G:{\mathbb{R}}^{n+m+1}\to {\mathbb{R}}^n$ the jump function, and $h:{\mathbb{R}}^{n+m+1}\to {\mathbb{R}}^p$ the output function, while $\mathcal{C}\subset {\mathbb{R}}^{n+m+1}$ represents the flow set and $\mathcal{D}\subset {\mathbb{R}}^{n+m+1}$ is the jump set. When executing an ode simulation, a jump condition trigger is permanently verified and, based on the selected configuration, it allows prioritizing the flow logic, the jump logic, or a stochastic behavior which includes randomly selecting any of them. This jump condition will also be needed for hybrid system interconnections.

We propose a separate class in the toolbox, called HybridSystem, which inherits the previously described class System, includes the ode event-based mechanism from HyEQ Toolbox [16], and is extended to support time-varying differential equation systems and exogenous input signals. Besides the base interface from System, it also provides methods for functions G, $\mathcal{C}$, and $\mathcal{D}$. It also provides a wrapper function to promote any System object to the type HybridSystem, by adding dummy G, $\mathcal{C}$, and $\mathcal{D}$ methods, in order to be compatible for use in hybrid system interconnections. The flexibility added by this class in the toolbox allows model-in-the-Loop simulations using physical processes with hybrid dynamics, such as switching systems, i.e., electrical machines and power converters, or simulations of the closed loop control system, seen as hybrid system through the interconnection of a continuous-time process and a discrete-time controller, allowing the user to assess several performance analysis steps.

2.3. System Interconnections

After defining individual or atomic systems as in previous sections, the necessity for composing system interconnections readily appears. The classical interconnection operations are the series, parallel, lower, and upper linear-fractional transformations (LLFT and ULFT). Moreover, two separate approaches have been considered, i.e., to interconnect general-purpose nonlinear systems modeled by the class System and hybrid systems modeled by the class HybridSystem separately. The first case is useful for linearization near an operating point, studying its system theoretical properties, and designing control techniques, while the latter becomes useful in a model-in-the-Loop simulation context and for closed loop system property analysis. All provided system interconnections are implemented in classes which inherit the base class System.

The software classes presented in this section extend the series, parallel, feedback, and lft functions from MATLAB for nonlinear and hybrid systems, based on the interfaces from Equations (1) and (4). For hybrid system interconnections, the continuous and discrete dynamics sets $\mathcal{C}$ and $\mathcal{D}$, respectively, are obtained using union and intersection set operations.

Moreover, the next discrete state for each subsystem is triggered by its own logic, predefined in the jump function G and only when necessary; otherwise, it remains unchanged. For specifying this next discrete state ${\mathbf{x}}^{+}$ logic, as in the interface from Equation (4c), we will use the notation IF(CONDITION, THEN, ELSE), where CONDITION will be true when the point in the state-space is in the jump set, i.e., $\left(\mathbf{x},\mathbf{u},t\right)\in \mathcal{D}$ or $\left(\mathbf{x},\mathbf{u},t\right)\notin \mathcal{C}$; THEN gives the next state if a jump needs to be performed; and ELSE gives the next discrete state otherwise.

The state, output, and hybrid domain equations for nonlinear and hybrid system series connection, with the notations used in Figure 2, upper row, implemented in classes SeriesConnectionSystem and HybridSeriesConnectionSystem, are as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& = \left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\\ F_2\left({\mathbf{x}}_2,h_1\left({\mathbf{x}}_1,\mathbf{u},t\right),t\right)\end{array}\right];\hfill \\ \hfill \mathbf{y}& = h_2\left({\mathbf{x}}_2,h_1\left({\mathbf{x}}_1,\mathbf{u},t\right),t\right).\hfill \end{array}{\displaystyle |}\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& = \left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\\ F_2\left({\mathbf{x}}_2,h_1\left({\mathbf{x}}_1,\mathbf{u},t\right),t\right)\end{array}\right];\hfill \\ \hfill \left[\begin{array}{c}{\mathbf{x}}_1^{+}\\ {\mathbf{x}}_2^{+}\end{array}\right]& = \left[\begin{array}{c}\mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{1}},G_1\left({\mathbf{x}}_1,\mathbf{u},t\right),{\mathbf{x}}_1\\ \mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{2}},G_2\left({\mathbf{x}}_2,h_1\left({\mathbf{x}}_1,\mathbf{u},t\right),t\right),{\mathbf{x}}_2\end{array}\right];\hfill \\ \hfill \mathcal{C}\left(\mathbf{x},\mathbf{u},t\right)& = {\mathcal{C}}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\cap {\mathcal{C}}_2\left({\mathbf{x}}_2,h_1\left({\mathbf{x}}_1,\mathbf{u},t\right),t\right);\hfill \\ \hfill \mathcal{D}\left(\mathbf{x},\mathbf{u},t\right)& = {\mathcal{D}}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\cup {\mathcal{D}}_2\left({\mathbf{x}}_2,h_1\left({\mathbf{x}}_1,\mathbf{u},t\right),t\right);\hfill \\ \hfill \mathbf{y}& = h_2\left({\mathbf{x}}_2,h_1\left({\mathbf{x}}_1,\mathbf{u},t\right),t\right).\hfill \end{array}\end{array}$$

(5)

Given two initial subsystems Sys1 and Sys2 with dimensions $(m_1, n_1, p_1)$ and $(m_2, n_2, p_2)$, respectively, the resulting series connection system will have dimensions $(m=m_1, n=n_1+n_2, p=p_2)$.

The state, output, and hybrid domain equations for nonlinear and hybrid system parallel connection, with the notations used in Figure 2, bottom row, implemented in the classes named ParallelConnectionSystem and HybridParallelConnectionSystem, are as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& = \left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\\ F_2\left({\mathbf{x}}_2,\mathbf{u},t\right)\end{array}\right];\hfill \\ \hfill \mathbf{y}& = h_1\left({\mathbf{x}}_1,\mathbf{u},t\right)+h_2\left({\mathbf{x}}_2,\mathbf{u},t\right).\hfill \end{array}{\displaystyle |}\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& = \left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\\ F_2\left({\mathbf{x}}_2,\mathbf{u},t\right)\end{array}\right];\hfill \\ \hfill \left[\begin{array}{c}{\mathbf{x}}_1^{+}\\ {\mathbf{x}}_2^{+}\end{array}\right]& = \left[\begin{array}{c}\mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{1}},G_1\left({\mathbf{x}}_1,\mathbf{u},t\right),{\mathbf{x}}_1\\ \mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{2}},G_2\left({\mathbf{x}}_2,\mathbf{u},t\right),{\mathbf{x}}_2\end{array}\right];\hfill \\ \hfill \mathcal{C}\left(\mathbf{x},\mathbf{u},t\right)& ={\mathcal{C}}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\cap {\mathcal{C}}_2\left({\mathbf{x}}_2,\mathbf{u},t\right);\hfill \\ \hfill \mathcal{D}\left(\mathbf{x},\mathbf{u},t\right)& ={\mathcal{D}}_1\left({\mathbf{x}}_1,\mathbf{u},t\right)\cup {\mathcal{D}}_2\left({\mathbf{x}}_2,\mathbf{u},t\right);\hfill \\ \hfill \mathbf{y}& = h_1\left({\mathbf{x}}_1,\mathbf{u},t\right)+h_2\left({\mathbf{x}}_2,\mathbf{u},t\right).\hfill \end{array}\end{array}$$

(6)

Given two initial subsystems Sys1 and Sys2 with dimensions $(m_1, n_1, p_1)$ and $(m_2, n_2,$$p_2)$, respectively, the resulting parallel connection system will have dimensions $(m=m_1=m_2, n=n_1+n_2, p=p_1=p_2)$.

The state, output, and hybrid domain equations for nonlinear and hybrid system lower linear fractional trasformation (LLFT) connection, with the notations used in Figure 3, upper row, implemented in classes LLFTConnectionSystem and HybridLLFTConnectionSystem, are as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& = \left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{LLFT}},t\right)\\ F_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{LLFT}},t\right)\end{array}\right];\hfill \\ \hfill \mathbf{y}& = \left[\begin{array}{c}{h}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{LLFT}},t\right)\\ h_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{LLFT}},t\right)\end{array}\right].\hfill \end{array}{\displaystyle |}\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& = \left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{LLFT}},t\right)\\ F_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{LLFT}},t\right)\end{array}\right];\hfill \\ \hfill \left[\begin{array}{c}{\mathbf{x}}_1^{+}\\ {\mathbf{x}}_2^{+}\end{array}\right]& =\left[\begin{array}{c}\mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{1}},G_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{LLFT}},t\right),{\mathbf{x}}_1\\ \mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{2}},G_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{LLFT}},t\right),{\mathbf{x}}_2\end{array}\right];\hfill \\ \hfill \mathcal{C}\left(\mathbf{x},\mathbf{u},t\right)& = {\mathcal{C}}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{LLFT}},t\right)\cap {\mathcal{C}}_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{LLFT}},t\right);\hfill \\ \hfill \mathcal{D}\left(\mathbf{x},\mathbf{u},t\right)& = {\mathcal{D}}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{LLFT}},t\right)\cup {\mathcal{D}}_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{LLFT}},t\right);\hfill \\ \hfill \mathbf{y}& = \left[\begin{array}{c}{h}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{LLFT}},t\right)\\ h_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{LLFT}},t\right)\end{array}\right],\hfill \end{array}\end{array}$$

(7)

with the predefined notations:  

$$\mathbf{u}=\left[\begin{array}{c}{\mathbf{u}}_1\\ {\mathbf{u}}_2\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_{11}\\ {\mathbf{u}}_{12}^{ref}\\ {\mathbf{u}}_{21}^{ref}\\ {\mathbf{u}}_{22}\end{array}\right],{\mathbf{u}}_1^{\mathtt{LLFT}}=\left[\begin{array}{c}{\mathbf{u}}_{11}\\ {\mathbf{u}}_{12}^{ref}+{\mathbf{y}}_{21}^{-}\end{array}\right]\equiv \left[\begin{array}{c}{\mathbf{u}}_{11}\\ {\mathbf{u}}_{12}\end{array}\right],{\mathbf{u}}_2^{\mathtt{LLFT}}=\left[\begin{array}{c}{\mathbf{u}}_{21}^{ref}+{\mathbf{y}}_{12}^{-}\\ {\mathbf{u}}_{22}\end{array}\right]\equiv \left[\begin{array}{c}{\mathbf{u}}_{21}\\ {\mathbf{u}}_{22}\end{array}\right].$$

(8)

The common convention in the literature is to consider the last NCON values from the input vector ${\mathbf{u}}_1$, i.e., ${\mathbf{u}}_{12}$, as control input signals, while the last NMEAS values from the output vector ${\mathbf{y}}_1$, i.e., ${\mathbf{y}}_{12}$, as measurements signals. Only the vector $\mathbf{u}$ will be an exogenous signal, as the feedback components ${\mathbf{y}}_{12}^{-}$ and ${\mathbf{y}}_{21}^{-}$ are local and private feedback components computed implicitly at the previous time step, dictated by the selected ode solver. The exogenous signals ${\mathbf{u}}_{11}$ and ${\mathbf{u}}_{22}$ are seen as disturbance signals, while the signals ${\mathbf{u}}_{12}^{ref}$ and ${\mathbf{u}}_{21}^{ref}$ are seen as reference signals.

Given two initial subsystems Sys1 and Sys2 with dimensions $\left(m_1, n_1, p_1\right)$ and $(m_2, n_2,$$p_2)$, respectively, the resulting LLFT connection system will have dimensions $(m=m_1+m_2, n=n_1+n_2, p=p_1+p_2)$. The subsystem Sys1 is usually seen as the controlled plant, while Sys2 is seen as the controller. In order to assure compatibility between the two, several assertions must be made: NMEAS = length(${\mathbf{y}}_{12}$) = length(${\mathbf{u}}_{21}$) and NCON = length(${\mathbf{y}}_{21}$) = length(${\mathbf{u}}_{12}$).

The state, output, and hybrid domain equations for nonlinear and hybrid system upper linear fractional trasformation (ULFT) connection, with the notations used in Figure 3, bottom row, implemented in classes ULFTConnectionSystem and HybridULFTConnectionSystem, are as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& =\left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{ULFT}},t\right)\\ F_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{ULFT}},t\right)\end{array}\right];\hfill \\ \hfill \mathbf{y}& =\left[\begin{array}{c}{h}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{ULFT}},t\right)\\ h_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{ULFT}},t\right)\end{array}\right].\hfill \end{array}{\displaystyle |}\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\end{array}\right]& =\left[\begin{array}{c}{F}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{ULFT}},t\right)\\ F_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{ULFT}},t\right)\end{array}\right];\hfill \\ \hfill \left[\begin{array}{c}{\mathbf{x}}_1^{+}\\ {\mathbf{x}}_2^{+}\end{array}\right]& =\left[\begin{array}{c}\mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{1}},G_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{ULFT}},t\right),{\mathbf{x}}_1\\ \mathtt{if}{\left(\mathtt{jump}\right.}_{\mathtt{2}},G_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{ULFT}},t\right),{\mathbf{x}}_2\end{array}\right];\hfill \\ \hfill \mathcal{C}\left(\mathbf{x},\mathbf{u},t\right)& ={\mathcal{C}}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{ULFT}},t\right)\cap {\mathcal{C}}_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{ULFT}},t\right);\hfill \\ \hfill \mathcal{D}\left(\mathbf{x},\mathbf{u},t\right)& ={\mathcal{D}}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{ULFT}},t\right)\cup {\mathcal{D}}_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{ULFT}},t\right);\hfill \\ \hfill \mathbf{y}& =\left[\begin{array}{c}{h}_1\left({\mathbf{x}}_1,{\mathbf{u}}_1^{\mathtt{ULFT}},t\right)\\ h_2\left({\mathbf{x}}_2,{\mathbf{u}}_2^{\mathtt{ULFT}},t\right)\end{array}\right],\hfill \end{array}\end{array}$$

(9)

with the predefined notations:

$$\mathbf{u}={\mathbf{u}}_{12},{\mathbf{u}}_1^{\mathtt{ULFT}}=\left[\begin{array}{c}{\mathbf{y}}_2^{-}\\ {\mathbf{u}}_{12}\end{array}\right],{\mathbf{u}}_2^{\mathtt{ULFT}}={\mathbf{y}}_{11}^{-}.$$

(10)

The common convention in the literature is to consider the first NU values from the input vector of the plant subsystem, i.e., ${\mathbf{u}}_{11}$, as input uncertainty signals, while the first NY values from the output vector of the plant, i.e., ${\mathbf{y}}_{11}$, as output uncertainty signals. Only the vector $\mathbf{u}\equiv {\mathbf{u}}_{12}$ will be an exogenous signal, as the feedback components ${\mathbf{y}}_{11}^{-}$ and ${\mathbf{y}}_2^{-}$ are local and private feedback components computed implicitly at the previous time step, dictated by the selected ode solver. The exogenous signal ${\mathbf{u}}_{12}$ is seen as set of performance and control signals for the plant, without any reference signals recalled explicitly compared to the LLFT case.

Given two initial subsystems Sys1 and Sys2 with dimensions $\left(m_1, n_1, p_1\right)$ and $(m_2, n_2,$$p_2)$, respectively, the resulting ULFT connection system will have dimensions $(m=m_1+m_2, n=n_1+n_2, p=p_1+p_2)$. The subsystem Sys1 is usually seen as the augmented controlled plant, while Sys2 is seen as the unstructured uncertainty block. In order to assure compatibility between the two, several assertions must be made: NY = length(${\mathbf{y}}_{11}$) = length(${\mathbf{u}}_2$) and NU = length(${\mathbf{y}}_2$) = length(${\mathbf{u}}_{11}$).

2.4. Automatic Equilibrium Point Computation

Given a nonlinear system as in Equation (1), which may also include interconnections of systems, an operating point is desired with some of the input, state, and output variables imposed, such as a water level in a tank $y_h\left(t\right)$ controlled through two pumps ${\mathbf{u}}_{flow}\left(t\right)$, one with variable flow and one with a fixed flow, or a mechanical transportation system having a desired velocity $y_{\omega }\left(t\right)$ with respect to input forces and loads $\mathbf{u}\left(t\right)$. As such, a mechanism to automatically compute a partially imposed equilibrium point for an entire family of uncertain plants, relative to one which is considered nominal at the design phase, is proposed in this paragraph.

Starting from the system definition with dimensions m, n, and p, consider the sets of indexes, denoted by $\mathcal{I}$, and prescribed values, denoted by $\mathcal{V}$, for the input, state, and output variables, respectively:

$$\{\begin{array}{}\mathrm{(11a)}& \hfill {\mathcal{L}}_{\mathbf{u}}:=& {\mathcal{I}}_{\mathbf{u}}\left\{i_1^u,i_2^u,\dots ,i_{{\overline{m}}_u}^u\right\},{\mathcal{V}}_{\mathbf{u}}:=\left\{\overline{\mathbf{u}}\left(i_1^u\right),\overline{\mathbf{u}}\left(i_2^u\right),\dots ,\overline{\mathbf{u}}\left(i_{{\overline{m}}_u}^u\right)\right\},0\le {\overline{m}}_u\le m;\hfill \mathrm{(11b)}& \hfill {\mathcal{L}}_{\mathbf{x}}:=& {\mathcal{I}}_{\mathbf{x}}\left\{i_1^x,i_2^x,\dots ,i_{{\overline{n}}_x}^x\right\},{\mathcal{V}}_{\mathbf{x}}:=\left\{\overline{\mathbf{x}}\left(i_1^x\right),\overline{\mathbf{x}}\left(i_2^x\right),\dots ,\overline{\mathbf{x}}\left(i_{{\overline{n}}_x}^x\right)\right\},0\le {\overline{n}}_x\le n;\hfill \mathrm{(11c)}& \hfill {\mathcal{L}}_{\mathbf{y}}:=& {\mathcal{I}}_{\mathbf{y}}\left\{i_1^y,i_2^y,\dots ,i_{{\overline{p}}_y}^y\right\},{\mathcal{V}}_{\mathbf{y}}:=\left\{\overline{\mathbf{y}}\left(i_1^y\right),\overline{\mathbf{y}}\left(i_2^y\right),\dots ,\overline{\mathbf{y}}\left(i_{{\overline{p}}_y}^y\right)\right\},0\le {\overline{p}}_y\le p,\hfill \end{array}$$

along with their complementary sets of values for the indexes, denoted by $\mathcal{UI}$, and the values, denoted $\mathcal{UV}$, to be computed through optimization by solving a system of equations:

$$\{\begin{array}{}\mathrm{(12a)}& {\mathcal{UI}}_{\mathbf{u}}:=\left\{i_1^u,i_2^u,\dots ,i_{{\tilde{m}}_u}^u\right\},{\mathcal{UV}}_{\mathbf{u}}:=\left\{\tilde{\mathbf{u}}\left(i_1^u\right),\tilde{\mathbf{u}}\left(i_2^u\right),\dots ,\tilde{\mathbf{u}}\left(i_{{\tilde{m}}_u}^u\right)\right\},0\le {\tilde{m}}_u\le m;\hfill \mathrm{(12b)}& {\mathcal{UI}}_{\mathbf{x}}:=\left\{i_1^x,i_2^x,\dots ,i_{{\tilde{n}}_x}^x\right\},{\mathcal{UV}}_{\mathbf{x}}:=\left\{\tilde{\mathbf{x}}\left(i_1^x\right),\tilde{\mathbf{x}}\left(i_2^x\right),\dots ,\tilde{\mathbf{x}}\left(i_{{\tilde{n}}_x}^x\right)\right\},0\le {\tilde{n}}_x\le n;\hfill \mathrm{(12c)}& {\mathcal{UI}}_{\mathbf{y}}:=\left\{i_1^y,i_2^y,\dots ,i_{{\tilde{p}}_y}^y\right\},{\mathcal{UV}}_{\mathbf{y}}:=\left\{\tilde{\mathbf{y}}\left(i_1^y\right),\tilde{\mathbf{y}}\left(i_2^y\right),\dots ,\tilde{\mathbf{y}}\left(i_{{\tilde{p}}_y}^y\right)\right\},0\le {\tilde{p}}_y\le p,\hfill \end{array}$$

with

$$\{\begin{array}{}\mathrm{(13a)}& {\overline{m}}_u+{\tilde{m}}_u=m,{\mathcal{I}}_u\cup {\mathcal{UI}}_u=\left\{1,2,\dots ,m\right\},{\mathcal{I}}_u\cap {\mathcal{UI}}_u=\varnothing ;\hfill \mathrm{(13b)}& {\overline{n}}_x+{\tilde{n}}_x=n,{\mathcal{I}}_x\cup {\mathcal{UI}}_x=\left\{1,2,\dots ,n\right\},{\mathcal{I}}_x\cap {\mathcal{UI}}_x=\varnothing ;\hfill \mathrm{(13c)}& {\overline{p}}_y+{\tilde{p}}_y=p,{\mathcal{I}}_y\cup {\mathcal{UI}}_y=\left\{1,2,\dots ,p\right\},{\mathcal{I}}_y\cap {\mathcal{UI}}_y=\varnothing .\hfill \end{array}$$

a set of permutation matrices $P_{\mathbf{u}}\in {\mathbb{R}}^{m\times m}$, $P_{\mathbf{x}}\in {\mathbb{R}}^{n\times n}$, $P_{\mathbf{y}}\in {\mathbb{R}}^{p\times p}$ are obtained after sorting the indexes such as the following system of vector-valued equations needs to be solved:

$$\{\begin{array}{}\mathrm{(14a)}& \mathbf{0}=F\left(P_{\mathbf{x}}·\left[\begin{array}{c}\overline{\mathbf{x}}\\ \tilde{\mathbf{x}}\end{array}\right],P_{\mathbf{u}}·\left[\begin{array}{c}\overline{\mathbf{u}}\\ \tilde{\mathbf{u}}\end{array}\right],\tilde{t}\right);\hfill \mathrm{(14b)}& P_{\mathbf{y}}·\left[\begin{array}{c}\overline{\mathbf{y}}\\ \tilde{\mathbf{y}}\end{array}\right]=h\left(P_{\mathbf{x}}·\left[\begin{array}{c}\overline{\mathbf{x}}\\ \tilde{\mathbf{x}}\end{array}\right],P_{\mathbf{u}}·\left[\begin{array}{c}\overline{\mathbf{u}}\\ \tilde{\mathbf{u}}\end{array}\right],\tilde{t}\right).\hfill \end{array}$$

The system from Equation (14) becomes equivalent to directly solving a system of equations of the form $\mathbf{0}=\mathcal{F}\left(\mathbf{z}\right)$ in the vector-valued unknown:

$$\mathbf{z}={\left[\begin{array}{ccccccc}{\tilde{\mathbf{x}}}^T& & {\tilde{\mathbf{u}}}^T& & {\tilde{\mathbf{y}}}^T& & \tilde{t}\end{array}\right]}^T\in {\mathbb{R}}^{{\tilde{n}}_x+{\tilde{m}}_u+{\tilde{p}}_y+1}.$$

(15)

If the dynamical system is time-invariant or if the required time value is known a priori, then the time variable can be removed from the solver or it can be imposed to a certain value $\overline{t}$ in the same manner as the for the other signals. Moreover, the method is flexible and allows imposing and solving only the subsystem (14a) if the output variables coincide with the states. The unknown variables from Equation (15) can be initialized to random values or a rough estimate for the entire family of uncertain plants can be obtained with the simulation to a step response of the nominal plant at the required amplitudes. All systems in the uncertain physical plant set will be found in the same mathematical vicinity. After solving the preferred algebraic system configuration from Equation (14), the desired equilibrium point $\left(\overline{\mathbf{u}},\overline{\mathbf{x}},\overline{\mathbf{y}},\overline{\mathbf{t}}\right)$ can be reconstructed using the inverse permutation matrices $P_{\mathbf{u}}^{-1}$, $P_{\mathbf{x}}^{-1}$, $P_{\mathbf{y}}^{-1}$ and the notations from Equations (11) and (12). The method findEqPoint from class System forms and solves the system (14) and computes the desired equilibrium point for its predefined dynamical system based on the specifications from Equations (11) and (12) given in the structure eqOpts.

After acquiring the desired equilibrium point, the system linearization can be easily deduced through numeric differentiation methods. The most straightforward method is to compute the first-order Jacobian matrices of the functions F and h with respect to the state and input signals $\mathbf{x}$ and $\mathbf{u}$, respectively:

$$A=\frac{\delta F}{\delta \mathbf{x}}{|}_{\left({\mathbf{x}}_0,{\mathbf{u}}_0,{\mathbf{t}}_0\right)};B=\frac{\delta F}{\delta \mathbf{u}}{|}_{\left({\mathbf{x}}_0,{\mathbf{u}}_0,{\mathbf{t}}_0\right)};C=\frac{\delta h}{\delta \mathbf{x}}{|}_{\left({\mathbf{x}}_0,{\mathbf{u}}_0,{\mathbf{t}}_0\right)};D=\frac{\delta h}{\delta \mathbf{u}}{|}_{\left({\mathbf{x}}_0,{\mathbf{u}}_0,{\mathbf{t}}_0\right)}.$$

(16)

The method linearize(${\mathbf{x}}_0$,${\mathbf{u}}_0$,$t_0$) from class System of Section 1 computes the matrices from Equation (16) with a first-order derivative approximation and also the output equilibrium value ${\mathbf{y}}_0=h\left({\mathbf{x}}_0,{\mathbf{u}}_0,t_0\right)$:  

$$A_{\overline{1,n},i}\approx \frac{F\left({\mathbf{x}}_0+\Delta {\mathbf{x}}_0^i,{\mathbf{u}}_0,t_0\right)-F\left({\mathbf{x}}_0,{\mathbf{u}}_0,t_0\right)}{\Delta x};B_{\overline{1,m},j}\approx \frac{F\left({\mathbf{x}}_0,{\mathbf{u}}_0+\Delta {\mathbf{u}}_0^j,t_0\right)-F\left({\mathbf{x}}_0,{\mathbf{u}}_0,t_0\right)}{\Delta u},$$

(17)

following that the output matrices C and D to be computed in a similar manner by replacing F with h in the above formulas. The shorthand notations are $\Delta {\mathbf{x}}_0^i={\left[\begin{array}{ccccccc}0& \cdots & 0& \Delta x& 0& \cdots & 0\end{array}\right]}^T$ and $\Delta {\mathbf{u}}_0^j={\left[\begin{array}{ccccccc}0& \cdots & 0& \Delta u& 0& \cdots & 0\end{array}\right]}^T$ for the disturbance vectors corresponding to the state with the index $i\in \overline{1,n}$ or input with index $j\in \overline{1,m}$, respectively, while the optimal [18] unit perturbations are, using double precision, $\Delta x=\mathtt{tol}·\left(1+\left|\right|{\mathbf{x}}_0\left|\right|\right)$ and $\Delta u=\mathtt{tol}·\left(1+\left|\right|{\mathbf{u}}_0\left|\right|\right)$, tol = ${10}^{-5}$. Obviously, when linearizing the system, the static amplification of the initial nonlinear system is not accounted in the procedure, but will not be relevant in the actual control design process and implementation due to the consideration of only Lipschitz function-based systems and, as such, it will be correctly compensated. The correct simulation of the linearized system near the operating point is done using the class LTIEqSystem.

2.5. Automatic Least Conservative Uncertainty Bound Computation

Figure 4 encompasses an overview of the toolbox classes described in Section 2.5 and Section 2.6, along with showing their relationship with the classes from the previous sections. Based on any desirable combination of the above system classes, i.e., System, LTISystem, LTIEqSystem, HybridSystem, and their interconnections, we propose a new functionality to aid in uncertain system modeling, ready for use in augmenting the plant for $\mu$ synthesis, implemented in the classes UncertainPlantFactory, seen in Figure 1, along with UncertaintyBoundOptimizationProblem, seen in Figure 4.

The common uncertainty model types considered in practice are gathered in

Table 1. Commonly used classes of perturbation and uncertainty models for multiple-input and multiple-output (MIMO) systems, implemented in class UncertainPlantFactory.

Uncertainty Type	Definition	Implementation
Additive	$G\left(s\right)=G_n\left(s\right)+\Delta \left(s\right)$	$\Delta \left(s\right)=G\left(s\right)-G_n\left(s\right)$
Inverse additive	${\left(G\left(s\right)\right)}^{-1}={\left(G_n\left(s\right)\right)}^{-1}+\Delta \left(s\right)$	$\Delta \left(s\right)={\left(G\left(s\right)\right)}^{-1}-{\left(G_n\left(s\right)\right)}^{-1}$
Input multiplicative	$G\left(s\right)=G_n\left(s\right)\left[I+\Delta \left(s\right)\right]$	$\Delta \left(s\right)={\left(G_n\left(s\right)\right)}^{-1}\left(G\left(s\right)-G_n\left(s\right)\right)$
Output multiplicative	$G\left(s\right)= \left[I+\Delta \left(s\right)\right]G_n\left(s\right)$	$\Delta \left(s\right)= \left(G\left(s\right)-G_n\left(s\right)\right){\left(G_n\left(s\right)\right)}^{-1}$
Inverse input multiplicative	${\left(G\left(s\right)\right)}^{-1}= \left[I+\Delta \left(s\right)\right]{\left(G_n\left(s\right)\right)}^{-1}$	$\Delta \left(s\right)={\left(G\left(s\right)\right)}^{-1}\left(G_n\left(s\right)\right)-I$
Inverse output multiplicative	${\left(G\left(s\right)\right)}^{-1}= {\left(G_n\left(s\right)\right)}^{-1}\left[I+\Delta \left(s\right)\right]$	$\Delta \left(s\right)= \left(G_n\left(s\right)\right){\left(G\left(s\right)\right)}^{-1}-I$
Left coprime factor	$G\left(s\right)= {\left(\tilde{M}+{\Delta }_{\tilde{M}}\right)}^{-1}\left(\tilde{N}+{\Delta }_{\tilde{N}}\right)$	${\Delta }_{\tilde{M}}= \left(\tilde{N}+{\Delta }_{\tilde{N}}\right){\left(G\left(s\right)\right)}^{-1}-\tilde{M},\mathbf{1}-\mathbf{DOF}$
Right coprime factor	$G\left(s\right)= \left(N+{\Delta }_N\right){\left(M+{\Delta }_M\right)}^{-1}$	${\Delta }_M= {\left(G\left(s\right)\right)}^{-1}\left(N+{\Delta }_N\right)-M,\mathbf{1}-\mathbf{DOF}$

. Besides the definition of the uncertain plant $G\left(s\right)$ starting from a nominal model $G_n\left(s\right)$ in relation to the uncertainty block $\Delta \left(s\right)$, the mathematical expression of $\Delta \left(s\right)$ is necessary to experimentally deduce its frequency response. Left and right coprime factor uncertainties are described by two blocks: ${\Delta }_M$ and ${\Delta }_N$, and one of them can be selected as a free term, i.e., one degree of freedom (1-DOF). The class UncertainPlantFactory provides an interface to define a nominal plant and a random plant from a prespecified set. Besides the methods getNominalPlant and getRandomPlant, both returning a System object, it implements each uncertainty type $\Delta \left(s\right)$ from Table 1, obtained through Monte Carlo simulation using the magnitude characteristic in the frequency domain.

The frequency domain relevant for the studied plant is defined in logarithmic scale as

$$\mathsf{\Omega }=\left\{\underline{\omega }={\omega }_1<{\omega }_2<\dots <{\omega }_{N-1}<{\omega }_N=\overline{\omega }\right\},$$

(18)

along with the magnitude characteristic values sampled at points from $\mathsf{\Omega }$. By convention, ${\left|\right|\Delta \left(j\omega \right)\left|\right|}_{\infty }\le 1,\forall \omega \ge 0$; thus, the worst-case uncertainty deduced experimentally through Monte Carlo simulation is written as $\Delta \left(s\right)W_{exp}\left(s\right)$. As an example for the additive uncertainty type, the uncertain plant family will be $G\left(s\right)=G_n\left(s\right)+\Delta \left(s\right)·W_{exp}\left(s\right)$, with ${\left|\right|\Delta \left(j\omega \right)\left|\right|}_{\infty }\le 1,\forall \omega \ge 0$ and the toolbox returns the worst-case experimental magnitude characteristic of $W_{exp}\left(s\right)$, sampled at points from $j\mathsf{\Omega }$ through the high-level method getUncertaintyModel-Boundary from class UncertainPlantFactory, which wraps over low-level methods such as getAdditiveUncBoundary, getInverseAdditiveUncBoundary etc.

Due to the fact that the sampled points from the previous paragraph cannot be directly accounted for in robust control synthesis and, moreover, they may not represent an actual transfer function, it appears the need to compute a least conservative low-order transfer function to model the desired uncertainty family. This problem has been solved by employing a global optimization algorithm, such as PSO, described in [19]. The PSO algorithm has been considered in favor of other global optimization algorithms, such as GA, due to its inherent structure of addressing semi-continuous functions, such as our strongly nonlinear semi-continuous function described in Equation (21) in the variable from Equation (19).

A particle $\mathbf{x}={\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_{n_4}\end{array}\right]}^T$ of the optimization problem is defined as the transfer function  

$$W_x\left(s\right)=\frac{k}{s^p}·\frac{\prod \left(\stackrel{\circ }{T}s+1\right)\prod \left(\frac{s^2}{{\stackrel{\circ }{\omega }}_n^2}+\frac{2\stackrel{\circ }{\zeta }}{{\stackrel{\circ }{\omega }}_n}s+1\right)}{\prod \left(\widehat{T}s+1\right)\prod \left(\frac{s^2}{{\widehat{\omega }}_n^2}+\frac{2\widehat{\zeta }}{{\widehat{\omega }}_n}s+1\right)}=\frac{x_1}{s^p}·\frac{{\displaystyle \prod _2^{n_1}\left(x_{k}s+1\right)\prod _{n_1+1}^{n_2}\left(\frac{s^2}{x_k^2}+\frac{2x_{k+1}}{x_k}s+1\right)}}{{\displaystyle \prod _{n_2+1}^{n_3}\left(x_{k}s+1\right)\prod _{n_3+1}^{n_4}\left(\frac{s^2}{x_k^2}+\frac{2x_{k+1}}{x_k}s+1\right)}},$$

(19)

with optimization variables

$$\begin{array}{c}k>0;\left\{{\stackrel{\circ }{T}}_i\in \left[\frac{1}{\overline{\omega }},\frac{1}{\underline{\omega }}\right]\right\},i\in \overline{2:n_1};\left\{{\widehat{T}}_i\in \left[\frac{1}{\overline{\omega }},\frac{1}{\underline{\omega }}\right]\right\},i\in \overline{n_2+1:n_3};\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill & \left\{\left({\stackrel{\circ }{\zeta }}_i\in \left(0,1\right),{\stackrel{\circ }{\omega }}_{n,i+1}\in \left[\underline{\omega },\overline{\omega }\right]\right)\right\},i\in \overline{n_1+1:2:n_2};\hfill \end{array}$$

(20b)

$$\begin{array}{cc}\hfill & \left\{\left({\widehat{\zeta }}_i\in \left(0,1\right),{\widehat{\omega }}_{n,i+1}\in \left[\underline{\omega },\overline{\omega }\right]\right)\right\},i\in \overline{n_3+1:2:n_4}.\hfill \end{array}$$

(20c)

The cost functional found to provide good results for most initial particle swarms is

$${\displaystyle J_{\mathsf{\Omega },W_{exp}}\left(W_x\right)=\sum _{k=1}^N\left||W_{exp}\left(j{\omega }_k\right){|}^{dB}-{\left|W_x\left(j{\omega }_k\right)\right|}^{dB}\right|·{\phi }_N\left(k\right)+\sum _{|W_x\left(j{\omega }_k\right)|<|W_{exp}\left(j{\omega }_k\right)|}\lambda ,}$$

(21)

where ${\phi }_N:\left\{1,2,\dots ,N\right\}\to {\mathbb{R}}^{+}$ from the first sum is a windowing function, based on the Gaussian window function, meant to amplify the penalization for low- and high-frequency components, while the second sum adds a high-cost term $\lambda ={10}^9$ for each frequency point $j{\omega }_k$ for which the candidate function $W_x$ is below the experimental reference uncertainty weight $W_{exp}$. The windowing function considered is defined by

$${\phi }_N\left(k\right)={\left(3-2·w_N\left(k-\frac{N+1}{2}\right)\right)}^2,$$

(22)

with the Gaussian function $w_N:\mathbb{R}\to {\mathbb{R}}^{+}$ and value $\alpha =2.5$:

$$w_N\left(x\right)=e^{-\frac{1}{2}{\left(\alpha \frac{x}{(N-1)/2}\right)}^2}=e^{-\frac{x^2}{2{\sigma }^2}}\iff log{w}_N\left(x\right)=-\frac{1}{2}{\left(\alpha \frac{x}{(N-1)/2}\right)}^2=-\frac{x^2}{2{\sigma }^2}.$$

(23)

The cost functional $J\left(W_x\right)$ is configurable, such that a different windowing function may be provided, or that the difference between $W_{exp}$ and $W_x$ in the integral may be considered in absolute magnitude values instead of decibels. An alternative formulation of this minimization problem would be to use a GA, where the optimization variable $\mathbf{x}=W_x\left(s\right)$ could mutate in order to obtain different transfer function structures: add or remove real poles and zeros or, also, complex conjugate pole and zero pairs.

This functionality is implemented in the class UncertaintyBoundOptimizationProblem, using the methods getTransferFunctionCandidate for Equation (19), computeCandidate-Fitness and fitnessIntegral for Equation (21) based on the magnitudes of $|W_{exp}\left(j\omega \right)|$ and $|W_x\left(j\omega \right)|$ for $\omega \in \mathsf{\Omega }$ from Equation (18) and optimize to use the PSO algorithm to compute the best candidate transfer function $W_{x,optim}\left(s\right)$. The plot function has been overloaded to facilitate seeing the fitness function values in real-time and, moreover, the Bode plot for the best candidate $W_x\left(s\right)$ compared to the experimental data $W_{exp}\left(s\right)$.

2.6. Robust Synthesis and Closed Loop Validation

Classical solutions to the ${\mathcal{H}}_2$/${\mathcal{H}}_{\infty }$ problems are presented in [20,21,22,23,24] and others. The control problem is typically formulated for the nominal plant using the generalized framework depicted in Figure 5a.

This generalized plant is obtained by augmenting the physical process model with a set of mathematical signals which aid the optimization procedure and has the following structure:

$$\begin{array}{c}\left(\begin{array}{c}\mathbf{z}\left(t\right)\\ \mathbf{y}\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{P}_{11}& P_{12}\\ P_{21}& P_{22}\end{array}\right)\left(\begin{array}{c}\mathbf{w}\left(t\right)\\ \mathbf{u}\left(t\right)\end{array}\right),\mathrm{with}\left(\begin{array}{cc}{P}_{11}& P_{12}\\ P_{21}& P_{22}\end{array}\right)=\left(\begin{array}{ccc}A& B_1& B_2\\ C_1& D_{11}& D_{12}\\ C_2& D_{21}& D_{22}\end{array}\right)\end{array},$$

(24)

where $\mathbf{w}\in {\mathbb{R}}^{n_w}$ is the exogenous input vector, $\mathbf{u}\in {\mathbb{R}}^{n_u}$ is the control input vector, $\mathbf{z}\in {\mathbb{R}}^{n_z}$ is the error (output, performance) vector, and $\mathbf{y}\in {\mathbb{R}}^{n_y}$ is the measurement vector. The closed loop system is given by the LLFT interconnection of P and K:

$$P_{\mathbf{z}\mathbf{w}}=\mathtt{LLFT}(P,K)=P_{11}+P_{12}K{\left(I-P_{22}K\right)}^{-1}{P}_{21}.$$

(25)

For the nominal plant, the target of the robust control problem is to minimize the ${\mathcal{H}}_{\infty }$ norm using a stabilizing controller K, which can be written as

$$\underset{K\mathrm{stab}.}{min}{∥P_{\mathbf{z}\mathbf{w}}∥}_{\infty }=\underset{K\mathrm{stab}.}{min}\underset{\omega \in {\mathbb{R}}^{+}}{sup}\overline{\sigma }\left(P_{\mathbf{z}\mathbf{w}}\left(j\omega \right)\right),$$

(26)

obtaining a (sub)optimal value $\gamma$ by iteration, which minimizes the effects of the input vector $\mathbf{w}\left(t\right)$ as seen through the performance output vector $\mathbf{z}\left(t\right)$.

However, this problem ensures only nominal stability and nominal performance. However, the plant is a model of a physical process, having uncertainties. There are two types of uncertainties: unstructured, which illustrates neglected and unmodelled dynamics and which are represented by a full block $\Delta \in {\mathbb{R}}^{m\times m}$, and parametric, which are represented by $\delta I$, where $\delta$ is the maximum bound of the variable parameter. In a mixed-scenario, the following set is considered:  

$$\Delta =\left\{\Delta =\mathrm{diag}\left({\delta }_1{I}_{n_1},\cdots ,{\delta }_s{I}_{n_s},{\Delta }_1,\cdots ,{\Delta }_f\right)|{\delta }_k\in \mathbb{R},{\Delta }_j\in {\mathbb{R}}^{m_j\times m_j},k=\overline{1,s},j=\overline{1,f}\right\}.$$

(27)

In Figure 5b, the closed loop system containing a LLFT connection between plant P and controller K and an ULFT connection between plant P and uncertainty block $\Delta$ is presented. In this case, the generalized plant contains one extra input vector, i.e., disturbance inputs $\mathbf{d}\in {\mathbb{R}}^{n_d}$, and one extra output vector, i.e., disturbance outputs $\mathbf{v}\in {\mathbb{R}}^{n_v}$, giving the following structure:  

$$P_{\Delta }\left(s\right)=\left(\begin{array}{ccc}{P}_{vd}\left(s\right)& P_{vw}\left(s\right)& P_{vu}\left(s\right)\\ P_{zd}\left(s\right)& P_{zw}\left(s\right)& P_{zu}\left(s\right)\\ P_{yd}\left(s\right)& P_{yw}\left(s\right)& P_{yu}\left(s\right)\end{array}\right)\iff P_{\Delta }:\left(\begin{array}{c}\dot{\mathbf{x}}\left(t\right)\\ \mathbf{v}\left(t\right)\\ \mathbf{z}\left(t\right)\\ \mathbf{y}\left(t\right)\end{array}\right)=\left(\begin{array}{cccc}A& B_d& B_w& B_u\\ C_v& D_{vd}& D_{vw}& D_{vu}\\ C_z& D_{zd}& D_{zw}& D_{zu}\\ C_y& D_{yd}& D_{yw}& D_{yu}\end{array}\right)\left(\begin{array}{c}\mathbf{x}\left(t\right)\\ \mathbf{d}\left(t\right)\\ \mathbf{w}\left(t\right)\\ \mathbf{u}\left(t\right)\end{array}\right).$$

(28)

A mathematical tool used for studying the robustness is the structured singular value, defined for a square matrix $M\in {\mathbb{C}}^{N\times N}$ with respect to the set $\Delta$ as

$${\displaystyle {\mu }_{\Delta }\left(M\right)=\frac{1}{{\displaystyle \underset{\Delta \in \Delta }{min}\left\{\overline{\sigma }(\Delta )\right|det(I-M\Delta )=0\}}},}$$

(29)

if there exists $\Delta \in \Delta$ such that the matrix $I-M\Delta$ is rank deficient; otherwise, it is 0. For the system presented in Figure 5b, the structured singular value of $\mathtt{LLFT}(P,K)$, according to $\Delta$, can be defined as

$${\mu }_{\Delta }\left(\mathtt{LLFT}(P,K)\left(s\right)\right)=\underset{\omega \in {\mathbb{R}}_{+}}{sup}{\mu }_{\Delta }\left(\mathtt{LLFT}(P,K)\left(j\omega \right)\right).$$

(30)

Besides the classical ${\mathcal{H}}_2/{\mathcal{H}}_{\infty }$ techniques, the $\mu$ synthesis framework manages to design a controller that meets the robust stability and robust performance specifications. The robust stability implies that a certain controller manages to stabilize all processes described by the upper linear fractional transformation between plant and uncertainty block, while the robust performance means that the controller is able to impose the desired closed loop performance in the worst-case scenario. Based on the main loop theorem, a controller K meets the robust stability and robust performance if and only if the structural singular value of the lower linear fractional transformation with respect to $\Delta$ is smaller than 1. Therefore, the minimization problem can be written as

$$\underset{K\mathrm{stab}.}{inf}\underset{\omega \in {\mathbb{R}}_{+}}{sup}{\mu }_{\Delta }\left(\mathtt{LLFT}(P,K)\left(j\omega \right)\right),$$

(31)

which is not a convex problem. Additionally, the structural singular values are difficult to be explicitly computed. In order to solve this problem, the following upper bound is used [25]:

$${\mu }_{\Delta }\left(\mathtt{LLFT}(P,K)\left(j\omega \right)\right)\le \underset{D\in \mathbf{D}}{inf}\overline{\sigma }(D·\mathtt{LLFT}(P,K)\left(j\omega \right)·D^{-1}),$$

(32)

where the set $\mathit{D}$ is defined in relation to the uncertainty set $\Delta$ as  

$$\mathit{D}=\left\{\mathrm{diag}\left(D_1,\cdots ,D_s,d_1{I}_{m_1},\cdots ,d_f{I}_{m_f}\right)|D_k=D_k^{\top }\in {\mathbb{R}}^{n_k\times n_k},d_j>0,k=\overline{1,s},j=\overline{1,f}\right\}.$$

(33)

Now, using this bound, the solution of the initial non-convex problem can be practically approximated by solving the following quasi-convex problem:

$$\underset{K\mathrm{stab}.}{inf}\underset{\omega \in {\mathbb{R}}_{+}}{sup}\underset{D\in \mathbf{D}}{inf}\overline{\sigma }\left(D\left(j\omega \right)·\mathtt{LLFT}(P,K)\left(j\omega \right)·{\left(D\left(j\omega \right)\right)}^{-1}\right).$$

(34)

Finally, if the system D is fixed, the problem (34) is nothing but a ${\mathcal{H}}_{\infty }$ control problem, in this case called the K step. Furthermore, for a fixed controller K, the D scale step can be obtained by solving a Parrot problem for a desired set of frequencies $\mathsf{\Omega }=\{{\omega }_1,\cdots ,{\omega }_N\}$ using a LMI and then obtain a minimum phase system after performing an identification step. Using these, an iterative algorithm, based on alternative D-K iterations, manages to solve the $\mu$ synthesis problem. This procedure starts with $D=I$ and successively applies a K step and a D scaling step until a stopping criterion is reached.

Of great use in the controller design phase are the sensitivity, complementary sensitivity, and control effort functions, respectively, defined by

$$\begin{array}{cc}\hfill S& :={\left(I+GK\right)}^{-1};\hfill \end{array}$$

(35a)

$$\begin{array}{cc}\hfill T& :=GK{\left(I+GK\right)}^{-1};\hfill \end{array}$$

(35b)

$$\begin{array}{cc}\hfill R& :=K{\left(I+GK\right)}^{-1}=KS,\hfill \end{array}$$

(35c)

where G is the open loop model. The great advantage of considering this approach is that it allows sculpting the relevant loop functions to impose steady-state and transitory regime performances, which are specified for different frequency ranges, using adequately selected weighting functions. Besides the minimization from Equation (34), different constraints can be added to the optimization problem to obtain a compromise between S, $KS$, and T at various frequencies:

$$\begin{array}{cc}\hfill & \underset{K\mathrm{stab}.}{inf}\underset{\omega \in {\mathbb{R}}_{+}}{sup}\underset{D\in \mathbf{D}}{inf}\overline{\sigma }\left(D\left(j\omega \right)·\mathtt{LLFT}(P,K)\left(j\omega \right)·{\left(D\left(j\omega \right)\right)}^{-1}\right),\hfill \end{array}$$

(36a)

$$\begin{array}{cc}\hfill & \mathrm{such}\mathrm{that}{∥{\left(\begin{array}{ccccc}{W}_{S}S& & W_{T}T& & W_{KS}KS\end{array}\right)}^T∥}_{\infty }<1,\hfill \end{array}$$

(36b)

also known as the mixed-sensitivity closed loop shaping $\mu$ synthesis method. The approach considered in this first iteration of the toolbox uses closed loop shaping with $\mu$ synthesis for the controller. Using the work in [26] as a starting point, different frequency response specifications, directly correlated with desired time-response performances, can be imposed in the weighting functions. The sensitivity weighting function $W_S$ depends on four parameters as in

$$W_S\left(s\right)={\left(\frac{\frac{1}{M^{1/n}}s+{\omega }_B}{s+{\omega }_B{A}^{1/n}}\right)}^n,$$

(37)

where ${\omega }_B$ represents the imposed bandwidth of the system; M imposes the ${\mathcal{H}}_{\infty }$ norm of the sensitivity function, in order to limit the overshoot of the system; n imposes the slope of the sensitivity function for low frequencies; and A imposes the maximum allowed steady-state error.

On the other hand, the complementary sensitivity weighting function $W_T$ can be generally defined by the following structure, in a symmetrical manner compared to $W_S$:

$$W_T\left(s\right)={\left(\frac{s+{\omega }_{BT}}{A_T^{1/n}s+{\omega }_{BT}{M}_T^{1/n}}\right)}^n,$$

(38)

with ${\omega }_{BT}$ being the imposed bandwidth of the system; $M_T$ imposes the ${\mathcal{H}}_{\infty }$ norm of $T\left(s\right)$; n imposes the roll-off slope of the closed loop system, which should be directly coupled with sensor noise characteristics; and $A_T$ imposes the least required attenuation for high frequencies. In practice, the complementary sensitivity bandwidth ${\omega }_{BT}$ can be adapted to the characteristics of the sensor in order to account for high-frequency noise.

Finally, the control effort weighting function with desired specifications $M_0:=\left|W_{KS}\left(0\right)\right|$, $M_{\infty }:=\left|W_{KS}(\infty )\right|$ and $|W_{KS}(j·{\omega }_d)|=M_d$, $M_0<M_d<M_{\infty }$ can be synthesized by the following formula:

$$W_{KS}\left(s\right)=\frac{M_{\infty }s+M_0{\omega }_d\sqrt{\frac{M_{\infty }^2-M_d^2}{M_d^2-M_0^2}}}{s+{\omega }_d\sqrt{\frac{M_{\infty }^2-M_d^2}{M_d^2-M_0^2}}}.$$

(39)

A higher-order counterpart can be generalized as for the previous cases, but was not found necessary for the proposed case studies and other tested benchmark plants.

Class RobustControlSynthesisProblem uses the nominal linearized plant around a required operating point using the system (14) and options (11)–(13), with a specified uncertainty type from Table 1, modeled through class UncertaintyBoundOptimization-Problem, allows imposing closed loop performance specifications with frequency weights (37)–(39), and synthesizes controller solutions for the problem (36) which cover the robust stability and robust performance problems. Additionally, the class also allows controller postprocessing, using order-reduction methods to compute easily implementable controllers. The aforementioned controller synthesis problem is illustrated in Figure 6, while the resulting LTI controller can be used in a MiL simulation context using the interface from class LTIEqSystem, encompassing system (3).

Class ClosedLoopControlProblem gathers all data, options, computations, and results from the previously presented individual blocks. It is used to define the end-to-end control problem by aggregating the uncertain plant family with its corresponding least conservative bound optimization; the robust control synthesis procedure and metadata; and allows Model-in-the-Loop simulations with the nonlinear, linearized, and hybrid system models, with automatic validation of imposed performances for the nonlinear plant in the frequency and time-domain alike.

3. Results

To showcase the ease of use and functionalities of the proposed toolbox, a set of case studies will be illustrated for DC-to-DC power converter circuits in a unified manner to encompass modeling, control synthesis and performance validation, considering the topologies buck, boost, and single-ended primary-inductor converter (SEPIC). DC-to-DC converters have been considered as a case study due to their ubiquity in various practical domains and applications, as presented in [27], ranging from renewable energy, hybrid and electric vehicles, controlled power sources, and many more. They can be seen as a good benchmark for control systems, due to their switching behavior, nonlinear dynamics, and different tried and tested control methods, such robust techniques in [28], Lyapunov methods in [29], passivity theory in [30], or sliding mode control as in [31].

This section will be split in a subsection which presents the converter mathematical models, a subsection with a suggested workflow for a general purpose control problem, followed by a subsection with numerical results and simulations for each of the studied converter topologies. For the SEPIC converter, having the most highly nonlinear behavior, thus being more difficult to control, we will illustrate and detail all plots generated by the toolbox, while for the buck and boost circuits, for brevity, we will show only the relevant figures and maintain the mathematical results and discussions.

3.1. Mathematical Modeling

The nonideal step-down (buck), step-up (boost), and single-ended primary-inductor converter (SEPIC) circuits are presented in Figure 7, where each component is described as follows.

• $S_1$: switching device, usually a transistor, and $S_2$: switching device, usually a diode or transistor;
• L, $L_1$, $L_2$: converter inductors;
• C, $C_{in}$, $C_1$, $C_2$: converter capacitors;
• R: (variable) output load resistance;
• E: external source voltage;
• $r_L$, $r_{L_1}$, $r_{L_2}$: resistances associated with the inductors;
• $r_C$, $r_{C_{in}}$, $r_{C_1}$, $r_{C_2}$: capacitors parasitic resistances;
• $r_{D{S}_1}$, $r_{D{S}_2}$: resistances associated with the ON state of the switching devices (usually drain source);
• $V_{F_1}$, $V_{F_2}$: constant voltage drops associated with the conducting phase of $S_1$ and $S_2$;
• $\mu \in \left[0,1\right]$: normalized duty cycle applied to $S_1$; complementary to the PWM signal applied to $S_2$.

Each converter has one control switching device $S_1$, while the other one, $S_2$, will be complementary to the former. Although, typically $S_2$ is a diode, it is preferable to use two encapsulated transistors for $S_1$ and $S_2$. When working in continuous conduction mode (CCM), all converters will have an ON state, corresponding to $S_1$ being on and $S_2$ off, along with an OFF state, for its complementary behavior. The corresponding LTI models, for a constant load resistance R, for the ON and OFF states will be presented using the following structure with the external voltage seen as disturbance input $u\left(t\right)=E\left(t\right)$, the voltage drops $V_{F_1}$ and $V_{F_2}$ as constant DC inputs, states from the inductor currents and capacitor voltages, and the load resistor voltage as measured output:

$$\left(\begin{array}{cc}A& \overline{B}\\ C& \overline{D}\end{array}\right)=\left(\begin{array}{ccc}A& B& B_V\\ C& D& D_V\end{array}\right),$$

(40)

where:  

$$\left\{\begin{array}{c}\dot{\mathbf{x}}=A_{ON}\mathbf{x}+B_{ON}E+B_{V,ON}\left[V_{F_1}{V}_{F_2}\right];\hfill \\ \mathbf{y}=C_{ON}\mathbf{x}+D_{ON}E+D_{V,ON}\left[V_{F_1}{V}_{F_2}\right],\hfill \end{array}\right.\left\{\begin{array}{c}\dot{\mathbf{x}}=A_{OFF}\mathbf{x}+B_{OFF}E+B_{V,OFF}\left[V_{F_1}{V}_{F_2}\right];\hfill \\ \mathbf{y}=C_{OFF}\mathbf{x}+D_{OFF}E+D_{V,OFF}\left[V_{F_1}{V}_{F_2}\right].\hfill \end{array}\right.$$

(41)

The control variable is the duty cycle of the switching devices $\mu \left(t\right) \in \left[0,1\right]$. Using a convex combination of the ON and OFF equation systems from Equation (41), an averaged state-space nonlinear model of the process is obtained close to the hybrid model’s behavior given a sufficiently high PWM frequency:

$$\dot{\mathbf{x}}\left(t\right)=\mu \left(t\right)·{\mathbf{x}}_{ON}\left(t\right)+(1-\mu \left(t\right))·{\mathbf{x}}_{OFF}\left(t\right)\equiv F\left(\mathbf{x}\left(t\right),{\left[E\left(t\right),R\left(t\right),\mu \left(t\right)\right]}^T,t\right).$$

(42)

As such, an affine nonlinear system, with respect to $\mu$, with the state function F as above can be implemented by inheriting the class System from the toolbox. The disturbances which affect the system are the voltage source E and variable output load R, which are stochastic in nature, along with uncertainties of its components due to manufacturing tolerances, relevant on inductors and capacitors. As the toolbox easily allows using the output capacitor voltage or output load voltage with minor modifications, we used the resistor voltage as measurement variable due to its corresponding practical control use cases. By inheriting the class UncertainPlantFactory, a set of tolerances can be imposed on all relevant circuit parameters and, also, an LTI uncertain set can be automatically computed with the provided mechanisms. The equilibrium point will have only the steady-state values of $\overline{E}$, $\overline{R}$, and ${\overline{u}}_R$ imposed, with the following structure:

$$\left(\overline{\mathbf{u}},\overline{\mathbf{x}},\overline{\mathbf{y}}\right)=\left(\left[\overline{E},\overline{R},\overline{\mu }\right],\left[{\overline{i}}_{L_1},{\overline{u}}_{C_1},\dots \right],{\overline{u}}_R\right).$$

(43)

After synthesizing a robust controller, model-in-the-loop simulations would be desired for the averaged state-space models and, also, for a hybrid model description of the converters. For the hybrid approach, the class UncertainPlantFactory is inherited again, this time with individual plants of type HybridSystem. For the general approach for CCM, based on the structure from Equation (4), the input vector is comprised of $\mathbf{u}={\left[E,R,\mu \right]}^T$, the state vector is extended to $\mathbf{x}={\left[{\mathbf{z}}^T,q,\tau \right]}^T$, with $\mathbf{z}={\left[i_{L_1},u_{C_1},\dots \right]}^T$ being the physical continuous states; $q\in \left\{0,1\right\}$ the discrete state number, i.e., ON and OFF; and $\tau \in \left[0,T_{PWM}\right)$ the time values for a single PWM period. After each PWM period completion, the auxiliary time state $\tau$ is reinitialized to 0. The model description becomes, for the state-space description

$$\{\begin{array}{}\mathrm{(44a)}& \left[\begin{array}{c}\dot{\mathbf{z}}\left(t\right)\\ \dot{q}\left(t\right)\\ \dot{\tau }\left(t\right)\end{array}\right]\hfill & =\left[\begin{array}{c}(1-q)·\left(A_{ON}\mathbf{z}\left(t\right)+{\overline{B}}_{ON}\mathbf{u}\left(t\right)\right)+q·\left(A_{OFF}\mathbf{z}\left(t\right)+{\overline{B}}_{OFF}\mathbf{u}\left(t\right)\right)\\ 0\\ 1\end{array}\right],\left(\mathbf{x},\mathbf{u},t\right)\in \mathcal{C};\hfill \mathrm{(44b)}& \left[\begin{array}{c}{\mathbf{z}}^{+}\left(t\right)\\ q^{+}\left(t\right)\\ {\tau }^{+}\left(t\right)\end{array}\right]\hfill & =\left[\begin{array}{c}\mathbf{z}\left(t\right)\\ \left\{\begin{array}{c}1,\mathtt{if}q==\mathtt{0}\\ 0,\mathtt{if}q==\mathtt{1}\end{array}\right.\\ \left\{\begin{array}{c}\tau ,\mathtt{if}q==\mathtt{0}\\ 0,\mathtt{if}q==\mathtt{1}\end{array}\right.\end{array}\right],\left(\mathbf{x},\mathbf{u},t\right)\in \mathcal{D},\hfill \end{array}$$

while the flow, jump, and output functions are

$$\{\begin{array}{}\mathrm{(45a)}& \hfill \mathcal{C}\left(\mathbf{x},\mathbf{u},t\right)& =\left\{\left(\left(q==0\right)\wedge \left(\tau \le \mu \left(t\right)·T_{PWM}\right)\right)\vee \left(\left(q==1\right)\wedge \left(\tau >\mu \left(t\right)·T_{PWM}\right)\right)\right\};\hfill \mathrm{(45b)}& \hfill \mathcal{D}\left(\mathbf{x},\mathbf{u},t\right)& =\left\{\left(\left(q==0\right)\wedge \left(\tau >\mu \left(t\right)·T_{PWM}\right)\right)\vee \left(\left(q==1\right)\wedge \left(\tau >T_{PWM}\right)\right)\right\};\hfill \mathrm{(45c)}& \hfill \mathbf{y}\left(t\right)& =h\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),t\right).\hfill \end{array}$$

In order to encompass DCM regimes, the number of discrete states must be extended with new LTI blocks obtained by adding the mathematical constraint of canceling the diode $S_2$ voltage and, as such, the corresponding current signal for that branch, along with more sophisticated jump functions, flow and jump sets. For brevity, we will not insist on these extensions, although an example described for the boost converter can be found in [32].

3.2. Toolbox Workflow

A suggested end-to-end workflow for the toolbox can be summarized in the following steps, all of which should be run from intermediary methods of an instance of class ClosedLoopControlProblem:

• inherit class System to define the nonlinear model of the process as in Equation (1) and Figure 1;
• define equilibrium point specifications as in (11)–(13);
• inherit class UncertainPlantFactory and overload method getRandomPlant;
• using the desired operating point, linearize a set of systems using (16) and experimentally determine an uncertainty model $W_{exp}$ based on Table 1;
• define uncertainty options, including transfer function structure for all particles, and execute the methods from class UncertaintyBoundOptimizationProblem in order to minimize the functional $J_{\mathsf{\Omega },W_{exp}}\left(W_x\right)$ from Equation (21), obtaining $W_{x,opt}$;
• run optimization to compute the uncertainty weight $W_{unc}\left(s\right)$ as in Table 1;
• define robust control specifications $W_S$, $W_T$, and $W_{KS}$ as in Equations (37)–(39);
• synthesize robust $\mu$ controller based on Figure 6;
• apply order-reducing methods on the resulting controller;
• validate frequency and time-response performance specifications using the nonlinear system at operating point through Model-in-the-Loop simulations;
• optionally, inherit the classes HybridSystem and UncertainPlantFactory, respectively, to validate time-response performance specifications using the corresponding hybrid plant model at operating point; for DC-to-DC converter control, the last step should be adapted for CCM or DCM operation.

3.3. Numerical Results

We will briefly present the obtained results for the three converters using the approaches established in Section 3.1 and Section 3.2.

3.3.1. SEPIC Converter

The SEPIC converter state-space model for the ON state of switch $S_1$, as structured in Equation (41), is

$$\left(\begin{array}{cccccccc}-\frac{1}{r_{C_{in}}{C}_{in}}& 0& 0& 0& 0& \frac{1}{r_{C_{in}}{C}_{in}}& 0& 0\\ \frac{1}{L_1}& -\frac{r_{C_{in}}+r_{L_1}+r_{D{S}_1}}{L_1}& 0& \frac{r_{D{S}_2}}{L_1}& 0& 0& -\frac{1}{L_1}& 0\\ 0& 0& 0& \frac{1}{C_1}& 0& 0& 0& 0\\ 0& -\frac{r_{D{S}_1}}{L_2}& -\frac{1}{L_2}& -\frac{r_{D{S}_1}+r_{C_1}+r_{L_2}}{L_2}& 0& 0& \frac{1}{L_2}& 0\\ 0& 0& 0& 0& -\frac{1}{(R+r_{C_2})C_2}& 0& 0& 0\\ 0& 0& 0& 0& \frac{R}{R+r_{C_2}}& 0& 0& 0\end{array}\right),$$

(46)

while for the OFF state of the switch $S_1$ is  

$$\left(\begin{array}{cccccccc}-\frac{1}{r_{C_{in}}{C}_{in}}& 0& 0& 0& 0& \frac{1}{r_{C_{in}}{C}_{in}}& 0& 0\\ \frac{1}{L_1}& -\frac{r_{aux}}{L_1}& -\frac{1}{L_1}& \frac{r_{D{S}_2}+r_{C_2}}{L_1}& -\frac{1}{L_1}& 0& 0& -\frac{1}{L_1}\\ 0& \frac{1}{C_1}& 0& 0& 0& 0& 0& 0\\ 0& \frac{r_{D{S}_2}+r_{C_2}}{L_2}& 0& -\frac{r_{D{S}_2}+r_{C_2}+r_{L_2}}{L_2}& \frac{1}{L_2}& 0& 0& \frac{1}{L_2}\\ 0& \frac{R}{(R+r_{C_2})C_2}& 0& -\frac{R}{(R+r_{C_2})C_2}& -\frac{1}{(R+r_{C_2})C_2}& 0& 0& 0\\ 0& \frac{r_{C}R}{R+r_C}& 0& -\frac{r_{C}R}{R+r_C}& \frac{R}{R+r_{C_2}}& 0& 0& 0\end{array}\right),$$

(47)

with the auxiliary notation $r_{aux}=r_{L_1}+r_{C_1}+r_{D{S}_1}+r_{C_2}+r_{C_{in}}$.

The nominal SEPIC converter parameters and their tolerances are presented in

Table 2. Single-ended primary-inductor converter (SEPIC) converter parameters, values, and corresponding tolerances.

Param.	Val.	Tol.	Param.	Val.	Tol.
$L_1$	$2.57$ [mH]	$\pm 20\%$	$L_2$	1.71[mH]	$\pm 20\%$
$r_{L_1}$	130 [m$\mathsf{\Omega }$]	$\pm 10\%$	$r_{L_2}$	110 [m$\mathsf{\Omega }$]	$\pm 10\%$
$r_{D{S}_1}$	$0.01$ [$\mathsf{\Omega }$]	$\pm 10\%$	$r_{D{S}_2}$	80 [m$\mathsf{\Omega }$]	$\pm 10\%$
$C_1$	$4.7$ [μF]	$\pm 20\%$	$C_2$	$3.57$ [μF]	$\pm 20\%$
$r_{C_1}$	270 [m$\mathsf{\Omega }$]	$\pm 10\%$	$r_{C_2}$	350 [m$\mathsf{\Omega }$]	$\pm 10\%$
$C_{in}$	$3.57$ [μF]	$\pm 20\%$	$r_{C_{in}}$	270 [m$\mathsf{\Omega }$]	$\pm 10\%$
$V_{F_1}$	$0.2$ [V]	$\pm 10\%$	$V_{F_2}$	$0.62$ [V]	$\pm 10\%$

.

The desired operating point specifications are output signal $y\left(t\right)\equiv u_R\left(t\right)$ is at 400 [V], with nominal voltage source and load inputs $u_1\left(t\right)$ = E = 300 [V], $u_2\left(t\right)$ = R = 80 $\left[\mathsf{\Omega }\right]$. The initial guesses for the state equilibrium values where $\tilde{\mathbf{x}}=[30,0.5,30,-0.5,30]$ and $u_3\left(t\right)$ = $\mu$ = $0.57$ for the duty cycle control input. After computation, the actual equilibrium point is $\left(\mathbf{u},\mathbf{x},\mathbf{y}\right) = \left(\left[300,80,0.5788\right],\left[300,6.8711,297.722,-5,400\right],400\right)$.

An input multiplicative uncertainty model, i.e., $G\left(s\right)=G_n\left(s\right)\left[1+\Delta \left(s\right)W_{unc}\left(s\right)\right]$, ${\left|\right|\Delta \left|\right|}_{\infty }\le 1$, as in Table 1 has been automatically computed from input $u_3\left(t\right)$ to output $y_1\left(t\right)$, with tolerances $\pm 10$ [V] and $\pm 5$$\left[\mathsf{\Omega }\right]$ for inputs $u_1\left(t\right)=E$ and $u_2\left(t\right)=R$, based on 1000 Monte Carlo simulations. The relevant frequencies vary in the interval $\left[\underline{\omega }={10}^{-2},\overline{\omega }={10}^8\right]$, with 300 equally distributed samples in log domain. A successful set of hyperparameters for the particle swarm optimization algorithm is comprised of a swarm size of 1000, initial swarm span of ${10}^4$, minimum neighbors fraction of $0.9$, and inertia range of $\left[0.1,1.1\right]$, for a transfer function structure as in Equation (19), with a complex pole pair and a complex zero pair, resulting in

$$W_{unc}\left(s\right)=\frac{0.67275\left(s^2+941.1s+2.222\times {10}^5\right)}{s^2+147.3s+5.422\times {10}^7}.$$

(48)

The linearized SEPIC plant family is comprised of fourth-order stable systems, in minimal form, with four zeros, three of which are of nonminimum phase. The nominal model is

$$G_n\left(s\right)=4.1368\frac{(s+8.003\mathrm{e}+5)(-s+2.304\mathrm{e}+4)(s^2-717.4s+5.145\mathrm{e}+7)}{(s^2+2673s+3.749\mathrm{e}+7)(s^2+1339s+6.493\mathrm{e}+7)}.$$

(49)

The entire uncertainty family has the same structure with poles, zeros, and equilibrium points in the vicinity of the nominal counterparts. The uncertain SEPIC plant family structure and behavior, along with the PSO conservative bound computation are illustrated in Figure 8. Figure 8-1 illustrates the pole-zero plot for the linearized uncertain SEPIC converter family, Figure 8-2 shows the best particle frequency-response fit $W_x\left(s\right)$ as in Equation (19) on the right y axis, and, also, the best functional fit on the left y axis, Figure 8-3 and Figure 8-4 show the frequency response of the plant $G\left(s\right)$ family and uncertainty family $\Delta \left(s\right)W_{unc}\left(s\right)$, respectively, while Figure 8-5 illustrates the system states and outputs for a $2\%$ step disturbance relative to the equilibrium input value ${\mu }_0=0.5788$.

For controller synthesis, the preferred loop shaping specifications were selected to highly penalize the control effort near the system resonance, as the obtained controllers would be difficult to implement in practice, needing high sampling frequencies, i.e., $f_e>20$ [kHz]. As such, for the sensitivity function: ${\omega }_B$= 200 [rad/s], $A={10}^{-2}$, $M=2$, $n=1$; for the complementary sensitivity function: ${\omega }_{BT}=2000$ [rad/s], $A_T={10}^{-4}$, $M_T=2$, $n=2$; and for the control effort $M_0=100$, $M_{\infty }={10}^5$, $|W_{KS}(j·200)|=250$, resulting in  

$$W_S\left(s\right)=\frac{0.5s+200}{s+2},W_T\left(s\right)=\frac{s^2+4000s+4\times {10}^6}{1\times {10}^{-4}{s}^2+56.57s+8\times {10}^6},W_{KS}\left(s\right)=\frac{{10}^{5}s+8.729\times {10}^6}{s+8.729\times {10}^4}.$$

(50)

From the $\mu$ synthesis procedure, a controller of order 21 is obtained. After order reduction, the smallest controller which manages to assure all imposed specifications for the plant family, with a peak value ${\mu }_{\Delta }\left(\mathtt{LLFT}(P,K)\right)\le 0.8361<1$, is given by the third-order system

$$K_{red}^{SEPIC}=\left(\begin{array}{cccc}-1.997& 3.056& 3.227& -0.5018\\ -3.057& -2197& -6118& -0.3838\\ 3.225& 6118& -1.016\times {10}^4& 0.4055\\ -0.5018& 0.3838& 0.4055& 0\end{array}\right).$$

(51)

The controller design phase, order reduction, and frequency response closed loop performance of the reduced-order one for the uncertain plant family are illustrated in Figure 9 and Figure 10. In this case, the control system has very large stability margins, with a phase margin of $\approx 82.1{[}^{\circ }]$ and gain margin in the interval $\left[19.4,20.3\right]$ [dB]. Additionally, as specified by the $n=2$ and $A_T={10}^{-4}$ parameters of the complementary sensitivity weighting function, $W_T$, the closed loop control system mitigates sensor noise signals with a considered spectrum starting from ${\omega }_{BT}>2000$ [rad/s], using an initial roll-off of $-40$ [dB/dec], followed by an attenuation of at least four orders of magnitude. In the actual MiL simulations, the attenuation does not stop at the prescribed value, as the system manages to maintain at least a $-20$ [dB/dec] roll-off.

From Figure 11, the closed loop bandwidth can be observed as ${\omega }_B>200$ [rad/s], equivalent to a rise time less than ≈5 [ms], a negligible steady-state error of ≈10${}^{-2}\times (y_{ss}-y_0)=0.2$ [V], where $y_{ss}$ represents the steady state value of the system, relative to the desired equilibrium value $y_0$. Moreover, the system has no overshoot, and it behaves like a first-order low-pass filter by design. The nonlinear MiL simulation options are $N=50$ random plants from the uncertainty set, solver is ode15i, due to the difficulty of simulating the closed loop plant otherwise (it is numerically unstable), with a step on the reference signal of $5\%$ from its initial equilibrium value of approximately 400 [V] and a simulation time of $0.05$ [s]. Almost the same conditions apply to the hybrid MiL simulation, with the solver switched to ode113, as ode15i is not supported here, and a PWM period selected randomly for each experiment from a nominal value $T_{PWM}=17.5$ [s] with a $\pm 20\%$ fluctuation from one simulation to another. A comparison of the nonlinear and hybrid cases is presented also in Figure 11, where it can be seen that the transitory regime is practically identical, but for the hybrid case, a slightly lower command signal is necessary in steady state. Due to working with $u_R\left(t\right)$ instead of $u_C\left(t\right)$ only for the output signal, the current ripple is propagated into the measured voltage.

3.3.2. Buck Converter

The buck converter state-space models for the ON and OFF states of switch $S_1$, respectively, as in Equation (41), are  

$$\left(\begin{array}{ccccc}-\frac{r_p+r_{D{S}_1}}{L}& -\frac{R}{(R+r_C)L}& \frac{1}{L}& -\frac{1}{L}& 0\\ \frac{R}{(R+r_C)C}& -\frac{1}{(R+r_C)C}& 0& 0& 0\\ \frac{r_{C}R}{R+r_C}& \frac{R}{R+r_C}& 0& 0& 0\end{array}\right);\left(\begin{array}{ccccc}-\frac{r_p-r_{D{S}_2}}{L}& -\frac{R}{(R+r_C)L}& 0& 0& -\frac{1}{L}\\ \frac{R}{(R+r_C)C}& -\frac{1}{(R+r_C)C}& 0& 0& 0\\ \frac{r_{C}R}{R+r_C}& \frac{R}{R+r_C}& 0& 0& 0\end{array}\right),$$

(52)

with the auxiliary notation $r_p=r_L+\frac{r_{C}R}{R+r_C}=r_L+r_C\left|\right|R$.

The nominal buck converter parameters and their tolerances are presented in

Table 3. Buck and boost converter parameters, values, and corresponding tolerances.

Param.	Val.	Tol.	Param.	Val.	Tol.
L	40 [μH]	$\pm 20\%$	$r_L$	10 [m$\mathsf{\Omega }$]	$\pm 10\%$
C	600 [μF]	$\pm 20\%$	$r_C$	$0.2$ [$\mathsf{\Omega }$]	$\pm 10\%$
$r_{D{S}_1}$	$0.01$ [$\mathsf{\Omega }$]	$\pm 10\%$	$r_{D{S}_2}$	$0.01$ [$\mathsf{\Omega }$]	$\pm 10\%$
$V_{F_1}$	$0.2$ [V]	$\pm 10\%$	$V_{F_2}$	$0.2$ [V]	$\pm 10\%$

.

The desired operating point specifications are as follows: the output signal $y\left(t\right)\equiv u_R\left(t\right)$ is at 5 [V], with nominal voltage source and load inputs $u_1\left(t\right)$ = E = 12 [V], $u_2\left(t\right)$ = R = 15 $\left[\mathsf{\Omega }\right]$. The initial guesses for the state equilibrium values where $\tilde{\mathbf{x}}=[1.25,5]$ and $u_3\left(t\right)=\mu =0.5$ for the duty cycle control input. After computation, the actual equilibrium point is $\left(\mathbf{u},\mathbf{x},\mathbf{y}\right)=\left(\left[12,15,0.4335\right],\left[0.333,4.999\right],5\right)$.

An input multiplicative uncertainty model, i.e., $G\left(s\right)=G_n\left(s\right)\left[1+\Delta \left(s\right)W_{unc}\left(s\right)\right]$, ${\left|\right|\Delta \left|\right|}_{\infty }\le 1$, as in Table 1, has been automatically computed from input $u_3\left(t\right)$ to output $y_1\left(t\right)$, with tolerances $\pm 1$ [V] and $\pm 1$$\left[\mathsf{\Omega }\right]$ for inputs $u_1\left(t\right)=E$ and $u_2\left(t\right)=R$, respectively, based on 1000 Monte Carlo simulations. The relevant frequencies vary in the interval $\left[\underline{\omega }={10}^1,\overline{\omega }={10}^7\right]$, with 200 equally distributed samples in log domain. A successful set of hyperparameters for the particle swarm optimization algorithm is comprised of a swarm size of 1000, initial swarm span of ${10}^4$, minimum neighbors fraction of $0.9$, and inertia range of $\left[0.1,1.1\right]$ for a transfer function structure as in Equation (19), with a real pole and a real zero, resulting in

$$W_{unc}\left(s\right)=\frac{0.51758\left(s+510.5\right)}{s+2906}.$$

(53)

The linearized buck plant family is comprised of second-order stable systems, in minimal form, with one zero. The nominal model is

$$G_n\left(s\right)=\frac{59178(s+8333)}{s^2+5261s+4.114\times {10}^7}.$$

(54)

The entire uncertainty family has the same structure with poles, zeros, and equilibrium points in the vicinity of the nominal counterparts. The uncertain buck plant family structure and behavior, along with the PSO conservative bound computation are illustrated in Figure 12.

For controller synthesis, the weighting function specifications were for the sensitivity function ${\omega }_B=1200$ [rad/s], $A={10}^{-4}$, $M=2$, $n=1$; for the complementary sensitivity function ${\omega }_{BT}=12\times {10}^3$ [rad/s], $A_T={10}^{-4}$, $M_T=2$, $n=2$; and for the control effort $M_0=0.1$, $M_{\infty }=100$, $|W_{KS}(j·1200)|=2$, resulting in  

$$W_S\left(s\right)=\frac{0.5s+1200}{s+0.12},W_T\left(s\right)=\frac{s^2+24000s+1.44\times {10}^8}{1\times {10}^{-4}{s}^2+339.4s+2.88\times {10}^8},W_{KS}\left(s\right)=\frac{100s+6006}{s+6.006\times {10}^4}.$$

(55)

As specified by the $n=2$ and $A_T={10}^{-4}$ parameters of the complementary sensitivity weighting function, $W_T$, the closed loop control system mitigates sensor noise signals with a considered spectrum starting from ${\omega }_{BT}>$ 12,000 [rad/s], using an initial roll-off of $-40$ [dB/dec], followed by an attenuation of at least four orders of magnitude. From the $\mu$ synthesis, a controller of order 17 is obtained. After order reduction, the smallest controller which manages to assure all imposed specifications for the plant family, with a peak ${\mu }_{\Delta }\left(\mathtt{LLFT}(P,K)\right)\le 0.97<1$, is

$$K_{red}^{Buck}=\left(\begin{array}{cccc}-0.12& -0.003479& -0.4751& 12.68\\ -0.0001768& -2.304& -1.241\times {10}^4& -0.1827\\ -0.4611& 1.241\times {10}^4& -4.522\times {10}^4& 25.09\\ 12.68& 0.1838& 25.09& 0\end{array}\right).$$

(56)

The nonlinear MiL simulation options are $N=50$ random plants from the uncertainty set, with the ode23t solver, with a step on the reference signal of $5\%$ from its initial equilibrium value of approximately 5 [V] and a simulation time of $0.02$ [s]. The simulation conditions for the hybrid MiL case are identical, with an additional PWM period selected randomly for each experiment from a nominal value $T_{PWM}=17.5$ [s] with a $\pm 20\%$ fluctuation from one simulation to another. A comparison of the nonlinear and hybrid cases is presented in Figure 13, where it can be seen that the transitory regime is practically identical. Due to working with $u_R\left(t\right)$ instead of $u_C\left(t\right)$ only for the output signal, the current ripple is propagated into the measured voltage.

3.3.3. Boost Converter

The boost converter state-space models for the ON and OFF states of switch $S_1$, respectively, as in Equation (41), are  

$$\left(\begin{array}{ccccc}-\frac{r_L+r_{D{S}_1}}{L}& 0& \frac{1}{L}& -\frac{1}{L}& 0\\ 0& -\frac{1}{(R+r_C)C}& 0& 0& 0\\ 0& \frac{R}{R+r_C}& 0& 0& 0\end{array}\right);\left(\begin{array}{ccccc}-\frac{r_{aux}}{L}& -\frac{R}{(R+r_C)L}& \frac{1}{L}& 0& -\frac{1}{L}\\ \frac{R}{(R+r_C)C}& -\frac{1}{(R+r_C)C}& 0& 0& 0\\ \frac{r_{C}R}{R+r_C}& \frac{R}{R+r_C}& 0& 0& 0\end{array}\right),$$

(57)

with the auxiliary notation $r_{aux}=r_L+r_{D{S}_2}+\frac{r_{C}R}{R+r_C}=r_L+r_{D{S}_2}+r_C\left|\right|R$. The nominal boost converter parameters and their tolerances are presented in Table 3, as they correspond with the buck converter parameters.

The desired operating point specifications are as follows: the output signal $y\left(t\right)\equiv u_R\left(t\right)$ is at 24 [V], with nominal voltage source and load inputs $u_1\left(t\right)$ = E = 12 [V], $u_2\left(t\right)$ = R = 15 $\left[\mathsf{\Omega }\right]$. The initial guesses for the state equilibrium values where $\tilde{\mathbf{x}}=[3,24]$ and $u_3\left(t\right)=\mu =0.5$ for the duty cycle control input. After computation, the actual equilibrium point is $\left(\mathbf{u},\mathbf{x},\mathbf{y}\right)=\left(\left[12,15,0.5179\right],\left[0.3189,23.999\right],24\right)$.

An input multiplicative uncertainty model, i.e., $G\left(s\right)=G_n\left(s\right)\left[1+\Delta \left(s\right)W_{unc}\left(s\right)\right]$, ${\left|\right|\Delta \left|\right|}_{\infty }\le 1$, as in Table 1, has been automatically computed from input $u_3\left(t\right)$ to output $y_1\left(t\right)$, with tolerances $\pm 1$ [V] and $\pm 1$$\left[\mathsf{\Omega }\right]$ for inputs $u_1\left(t\right)=E$ and $u_2\left(t\right)=R$, respectively, based on 1000 Monte Carlo simulations. The relevant frequencies vary in the interval $\left[\underline{\omega }={10}^1,\overline{\omega }={10}^7\right]$, with 200 equally distributed samples in log domain. A successful set of hyperparameters for the particle swarm optimization algorithm is comprised of a swarm size of 1500, initial swarm span of ${10}^4$, minimum neighbors fraction of $0.9$, and inertia range of $\left[0.1,1.1\right]$ for a transfer function structure as in Equation (19), with two real poles and two real zeros, resulting in

$$W_{unc}\left(s\right)=\frac{0.26592\left(s+512.5\right)\left(s+3.535\times {10}^4\right)}{\left(s+4016\right)\left(s+1.389\times {10}^4\right)}.$$

(58)

The linearized boost plant family is comprised of second-order stable systems, in minimal form, with two zeros, one being of nonminimum phase. The nominal model is

$$G_n\left(s\right)=\frac{0.65505(-s+8.551\times {10}^4)(s+8333)}{s^2+2988s+9.746\times {10}^6}.$$

(59)

The entire uncertainty family has the same structure with poles, zeros, and equilibrium points in the vicinity of the nominal counterparts. The uncertain boost plant family structure and behavior, along with the PSO least conservative second-order bound computation are illustrated in Figure 14.

For controller synthesis, the weighting function specifications were for the sensitivity function ${\omega }_B=650$ [rad/s], $A={10}^{-4}$, $M=2$, $n=1$; for the complementary sensitivity function ${\omega }_{BT}=3250$ [rad/s], $A_T={10}^{-4}$, $M_T=2$, $n=1$; and for the control effort $M_0=0.1$, $M_{\infty }=100$, $|W_{KS}(j·650)|=2$, resulting in

$$W_S\left(s\right)=\frac{0.5s+650}{s+0.065},W_T\left(s\right)=\frac{s+3250}{0.0001s+6500},W_{KS}\left(s\right)=\frac{100s+3253}{s+3.252\times {10}^4}.$$

(60)

An intrinsic limitation for the boost converter is that the bandwidth must not be imposed more than half of the value of the right-half plane non-minimum phase zero of the process, i.e., ${\omega }_B\le \frac{z}{2}\approx$ 43,000 [rad/s]. For this problem, this is a sufficiently high margin, as we also imposed a limitation for the command signal through $W_{KS}$. As specified by the $n=2$ and $A_T={10}^{-4}$ parameters of the complementary sensitivity weighting function, $W_T$, the closed loop control system mitigates sensor noise signals with a considered spectrum starting from ${\omega }_{BT}>3250$ [rad/s], using an initial roll-off of $-40$ [dB/dec], followed by an attenuation of at least four orders of magnitude. From the $\mu$ synthesis, a controller of order 17 is obtained. After order reduction, the smallest controller which retains all imposed specifications for the plant family, with a peak value ${\mu }_{\Delta }\left(\mathtt{LLFT}(P,K)\right)\le 0.9547<1$, is

$$K_{red}^{Boost}=\left(\begin{array}{cccc}-0.065& 0.8025& -0.002966& 5.007\\ 0.7116& -9.715\times {10}^4& 1.066\times {10}^4& -30.91\\ 0.002575& -1.065\times {10}^4& -1.422& -0.113\\ 5.007& -30.91& 0.1142& 0\end{array}\right).$$

(61)

With this regulator, the boost converter control system with parameters from Table 3 has very large stability margins, with phase margins between $\left[81,101\right]$${[}^{\circ }]$ and gain margins in the interval $\left[40,46\right]$ [dB]. The obtained sensitivity bandwidths vary between $\left[828,2510\right]$ [rad/s], all of them better than the prespecified value of 800.

Nonlinear MiL simulation options are $N=50$ random plants from the uncertainty set, with the ode15i solver, with a step on the reference signal of $5\%$ from its initial equilibrium value of approximately 5 [V] and a simulation time of $0.02$ [s]. The simulation conditions for the hybrid MiL case are almost identical, with the use of the ode113 solver instead, set to a relative tolerance of ${10}^{-8}$, and with an additional PWM period selected randomly for each experiment from a nominal value $T_{PWM}=17.5$ [s] with a $\pm 20\%$ fluctuation from one simulation to another. A comparison of the nonlinear and hybrid cases is presented in Figure 15, where it can be seen that the transitory regime is practically identical. Due to working with $u_R\left(t\right)$ instead of $u_C\left(t\right)$ only for the output signal, the current ripple is propagated into the measured voltage.

4. Conclusions and Future Work

There are several aspects to be discussed with regards to the presented toolbox and, additionally, including future work to be implemented in upcoming iterations.

The main reason for which the hybrid system framework is considered for the proposed toolbox is that the continuous-time plant needs to be regulated with a numerical controller. The hybrid system context may include LTI or nonlinear systems, switching systems (which are hybrid by nature), and singular systems, represented using differential-algebraic equations (DAEs). As such, for a continuous-time plant, a continuous-time controller is designed using the robust control framework. The controller is required to be numerically implementable, therefore two interfaces are necessary, which lead to a hybrid system. The numerical implementation must be easily obtained and automatically validated using rapid control prototyping (RCP) techniques. However, the properties of the closed loop system need to be reanalyzed after the discretization of the controller. Differential-algebraic equations represent a useful framework for modeling dynamical systems in engineering with a network-based structure of components. They are used in various industry fields such as mechanics (e.g., multiple-link mobile manipulator model), chemical engineering (modeling of chemical reactions), electrical engineering, cyber-physical systems, etc. All categories of processes taken into consideration in the hybrid framework are described using an approximate model that incorporates their relevant behavior. However, these types of systems have model and structure uncertainties. Moreover, the nonlinear systems that are linearized around an equilibrium point introduce such uncertainties as well. Therefore, to consider these uncertainties for the controller design, the robust control framework is mandatory. The first main applicability of RCP was to derive the necessary C/C++ source code with drivers for a given target microprocessor, and to simulate in reproducible conditions the behavior of a complex system. For the latter case, the most relevant simulation types, given in increasing order of complexity and closeness to reality, are Model-in-the-Loop (MiL), Software-in-the-Loop (SiL), and Hardware-in-the-Loop (HiL).

Many modeling software programs return circuits or mechanical systems already in DAE form, and it would be difficult, or sometimes impossible, to reformulate them in an ODE form without changing variables and losing their intended physical significance. In the context of the robust control framework used for DAEs, a significant work is the monograph in [33]. Although the robust control theory was well formulated in the previous years, this is still an open domain for research and publication, as surveyed and described in [34]. The difficulty of using methods which work correctly for multiple operating points is mitigated by using adaptive methods, such as gain scheduling for tracking problems. Furthermore, other relevant problems are automatic C/C++ code generation for controller implementation and the commutation between the prescribed operating points. Considering use cases for modeling, simulation, computer-aided design, and RCP, a relevant set of examples in the domain of power electronics, hybrid vehicles, and renewable energy systems is given in [35,36], where, although the main limitation is that analysis, control, and implementation aspects must be performed individually for the presented applications, they can be generally included under the same software framework, with the important exception of the robust control design and verification, along with the quantization sensitivity analysis in a unitary manner.

The main hindrance is that the majority of applications involve DAE systems of index greater than 1, meaning they necessitate more than one derivation step in order to formulate the problem, so there is a great need of tools that can also deal with high order DAEs, i.e., index 2 or 3 [37,38,39]. Another reason for the proposal of considering DAEs model is that they are also explicitly related to control issues, regarding both physical and operational constraints. Two illustrative examples are the case of improper systems, such as an ideal PID controller, and the elements that realize the decoupling in the MIMO systems case and when the plant has impulsive dynamics. The robust control techniques’ drawback is that the closed loop ${\mathcal{H}}_2/{\mathcal{H}}_{\infty }$ norm must be minimized using the prescribed weighting functions that penalize the exogenous outputs and these weighting functions need to be found ad hoc, which sometimes lead to some intermediary bad controllers and work overhead. Furthermore, after numerical implementation, the discrete controller loses a part of the imposed performances due to an inadequate sampling period or badly selected quantization levels. Although there are solutions in the literature, there is no unified approach to solve all these mentioned problems. As an extension to the framework and mindset given by the two previously mentioned RCP use cases, the current project proposed a highly automated toolbox for robust control design, which, in the current state-of-the-art is a highly iterative design process, when taking into consideration plant uncertainties, although the mathematical background for solving the optimization problems is well established. It proposes to eliminate design overhead when considering and modifying a specification set, manually redesigning the weighting functions, the optimization procedure, discretization of the regulators, quantization analysis, and closed loop analysis for the linearized and initial hybrid plant. Furthermore, in unison with high-performance numerical toolboxes, a justified report should automatically result after its use and explicitly state when unrealistic design specifications were considered.