Computational Techniques for Model Predictive Control of Large-Scale Systems with Continuous-Valued and Discrete-Valued Inputs

Abstract

We propose computational techniques for model predictive control of large-scale systems with both continuous-valued control inputs and discrete-valued control inputs, which are a class of hybrid systems. In the proposed method, we introduce the notion of virtual control inputs, which are obtained by relaxing discrete-valued control inputs to continuous variables. In online computation, first, we find continuous-valued control inputs and virtual control inputs minimizing a cost function. Next, using the obtained virtual control inputs, only discrete-valued control inputs at the current time are computed in each subsystem. In addition, we also discuss the effect of quantization errors. Finally, the effectiveness of the proposed method is shown by a numerical example. The proposed method enables us to reduce and decentralize the computation load.

1. Introduction

Control of large-scale systems is one of the fundamental problems in control theory and has been extensively studied so far (see, e.g., [1–7]). From the viewpoint of total optimization, it is necessary to regard a set of individual systems as one system. In recent years, it is important to consider energy management systems such as smart grid [8] and HEMS (Home Energy Management System) [9–11] as an application. As the other example, we give air-conditioning systems [12, 13] and so on. Thus, it is also important at the current stage to consider control of large-scale systems.

In large-scale systems, many kinds of actuators are included. For example, in air-conditioning systems, the output of air conditioners can be regarded as a continuous-valued control input, and ON/OFF switches of ceiling fans can be regarded as a discrete-valued control input. Thus, it is appropriate to consider large-scale systems with both continuous-valued control inputs and discrete-valued control inputs. However, to our knowledge, such a large-scale system has not been directly considered so far. Furthermore, such a large-scale system is a class of hybrid systems, and the finite-time optimal control problem is in general reduced to a mixed integer quadratic programming (MIQP) problem. For large-scale systems, the computation time for solving the MIQP problem is too long, and it is difficult to realize the model predictive control (MPC) method in which the MIQP problem is solved at each time.

We, thus, propose computational techniques for model predictive control of large-scale systems with both continuous-valued control inputs and discrete-valued control inputs. In the proposed method, we introduce the notion of virtual control inputs, which are obtained by relaxing discrete-valued control inputs to continuous variables. This term has been used in, for example, [12], but we define a different notion. By using the virtual control input, the MIQP problem is approximately rewritten as a quadratic programming (QP) problem, which can be relatively solved faster than the MIQP problem.

In online computation, first, continuous-valued control inputs and virtual control inputs minimizing a cost function are found in a centralized controller. Next, using the obtained virtual control inputs, only discrete-valued control inputs at the current time are computed in each decentralized controller. In other words, a sequence of discrete-valued control inputs is not computed because in MPC only the control input at the current time is applied to the plant. In a centralized controller, the QP problem is solved at each time. In each decentralized controller, the problem of finding the discrete-valued control input is solved at each time. This problem is reduced to an integer programming problem, or can be solved by the look-up table method in which a look-up table is generated off-line. Thus, the computation load is reduced and decentralized by the proposed method.

In addition, we also discuss the effect of quantization errors based on our previous result [14]. The effectiveness of the proposed method is shown by a numerical example. The proposed method provides us an easy-to-use method in control of large-scale systems.

This paper is organized as follows. In Section 2, large-scale systems studied here are explained, and an example of air-conditioning systems is shown. In Section 3, the optimal control problem is formulated. In Section 4, after the outline of the proposed method is explained, the details are explained. In Section 5, we discuss quantization errors. In Section 6, we show a numerical example. In Section 7, we conclude this paper.

Notation. Let ℛ denote the set of real numbers. Let 𝒵 denote the set of integers. Let 𝐼_𝑛 and 0_𝑚×𝑛 denote the 𝑛×𝑛 identity matrix and the 𝑚×𝑛 zero matrix, respectively. For simplicity, we sometimes use the symbol 0 instead of 0_𝑚×𝑛, and the symbol 𝐼 instead of 𝐼_𝑛. For a matrix 𝑀, let the matrix |𝑀| denote a matrix such that each element is given by an absolute value of each element of 𝑀.

2. Large-Scale Systems

Consider the discrete-time large-scale system consisting of 𝑠 subsystems given by Σ_𝑖∶𝑥_𝑖(𝑘+1)=^𝑠_𝑗=1𝐴_𝑖𝑗𝑥_𝑗(𝑘)+𝐵^𝑐_𝑖𝑢^𝑐_𝑖(𝑘)+𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘), (1) where 𝑖=1,2,…,𝑠, 𝑘=0,1,2,… is the discrete-time, and 𝑥_𝑖(𝑘)∈ℛ^𝑛𝑖 is the state in the subsystem 𝑖. The control inputs 𝑢^𝑐_𝑖(𝑘)∈𝒰^𝑐_𝑖⊆ℛ^𝑚𝑐𝑖 and 𝑢^𝑑_𝑖(𝑘)∈ℛ^𝑚𝑑𝑖 are the continuous-valued control input and the discrete-valued control input in the subsystem 𝑖, respectively. The set 𝒰^𝑐_𝑖 is given as a closed convex set. The discrete-valued control input 𝑢^𝑑_𝑖(𝑘) is given by 𝑢^𝑑_𝑖(𝑘)=𝑢^𝑑_𝑖,1(𝑘)𝑢^𝑑_𝑖,2(𝑘)⋮𝑢^𝑑_{𝑖,𝑚𝑑𝑖}(𝑘),𝑢^𝑑_𝑖,𝑙(𝑘)∈𝑢^𝑑_𝑖,𝑙,1,𝑢^𝑑_𝑖,𝑙,2,…,𝑢^𝑑_{𝑖,𝑙,𝑝𝑖,𝑙}=∶𝒰^𝑑_𝑖,𝑙⊆ℛ¹, (2) where 𝒰^𝑑_𝑖,𝑙 is the finite set expressing the candidates of the 𝑙th element of the discrete-valued control input in the subsystem 𝑖, and 𝑝_𝑖,𝑙 is the number of elements of 𝒰^𝑑_𝑖,𝑙. In system (1), the pair (𝐴_𝑖𝑖,𝐵^𝑐_𝑖𝐵^𝑑_𝑖) expresses the dynamics of the subsystem 𝑖, and the matrices 𝐴_𝑖𝑗,𝑖≠𝑗 express the effect of other subsystems 𝑗,𝑖≠𝑗.

We show an example of air-conditioning systems.

Example 1. Suppose that the dynamics of temperature in room 𝑖 are given by the following first-order plus time delay model: 𝑃_𝑖(𝑠)=𝐾_𝑖1+𝑇_𝑖𝑠𝑒^−𝐿𝑠, (3) where 𝐾_𝑖, 𝑇_𝑖, and 𝐿 are the static process gain, the time constant, and the time delay, respectively. For simplicity of discussion, the value of 𝐿 does not depend on room 𝑖. From (3), we can obtain the following state equation: 𝑞̇_𝑖(𝑡)1=−𝑇_𝑖𝑞_𝑖(𝑡)+𝐾_𝑖𝑇_𝑖𝑢_𝑖(),𝑦𝑡−𝐿_𝑖(𝑡)=𝑞_𝑖(𝑡), (4) where 𝑞_𝑖 is the state that implies the temperature in one room, and 𝑢_𝑖 is the continuous-valued control input that implies, for example, the output of a heat pump. Hereafter, the output equation 𝑦_𝑖(𝑡)=𝑞_𝑖(𝑡) is omitted.

Next, consider transforming (4) into a discrete-time system. For simplicity of discussion, suppose that the sampling period is given as 𝐿. Then, we can obtain the following discrete-time linear system: 𝑞_𝑖(𝑘+1)=𝑒^{−𝐿/𝑇𝑖}𝑞_𝑖(𝑘)+^𝐿₀𝑒^{−𝜏/𝑇𝑖}𝐾_𝑖𝑇_𝑖𝑑𝜏𝑢_𝑖(𝑘−1)=𝑒^{−𝐿/𝑇𝑖}𝑞_𝑖(𝑘)+𝐾_𝑖1−𝑒^{−𝐿/𝑇𝑖}𝑢_𝑖(𝑘−1). (5) By defining 𝑣_𝑖(𝑘)∶=𝑢_𝑖(𝑘−1), we can obtain 𝑞_𝑖()𝑣𝑘+1_𝑖()=𝑒𝑘+1^{−𝐿/𝑇𝑖}𝐾_𝑖1−𝑒^{−𝐿/𝑇𝑖}𝑞00^𝑖(𝑘)𝑣_𝑖(𝑘)+01𝑢_𝑖(𝑘). (6) Furthermore, consider adding discrete-valued control inputs. We consider three kinds of discrete-valued inputs as follows: (i) 𝑢^𝑤_𝑖∈{0,1}: window (close/open), (ii) 𝑢^ℎ_𝑖∈{0,1,2}: local heater (off/low/high), and (iii) 𝑢^𝑣_𝑖∈{0,1,2}: ventilation fan (off/low/high). Then, we can obtain the following temperature model in room 𝑖: 𝑞_𝑖()𝑣𝑘+1_𝑖(𝑘+1_{𝑥𝑖()𝑘+1})=𝑒^{−𝐿/𝑇𝑖}𝐾_𝑖1−𝑒^{−𝐿/𝑇𝑖}00_𝐴𝑖𝑞_𝑖(𝑘)𝑣_𝑖(𝑘)_{𝑥𝑖(𝑘)}+0_𝐵1_𝑐𝑖𝑢_𝑖(𝑘)_{𝑢𝑐𝑖(𝑘)}+𝑏^𝑤_𝑖𝑏^ℎ_𝑖𝑏^𝑣_𝑖000_𝐵𝑑𝑖𝑢^𝑤_𝑖(𝑘)𝑢^ℎ_𝑖(𝑘)𝑢^𝑣_𝑖(𝑘)_{𝑢𝑑𝑖(𝑘)}, (7) where 𝑏^𝑤_𝑖≤0, 𝑏^ℎ_𝑖≥0, and 𝑏^𝑣_𝑖≤0 are constants and are given in advance. We assume that outside temperature is equal to or less than room temperature.

Finally, we consider the effect of other rooms. Assume that airflow is expressed as a directed graph. In this example, suppose that airflow is given as in Figure 1. That is, we can obtain the directed graph in Figure 2. In addition, Figure 1 implies that room 5 has no window. That is, 𝑏^𝑤₅=0 holds. From Figure 2, the coupled term from room 𝑖 to room 𝑗 can be obtained as 𝑎_𝑗𝑖000_𝐴𝑗𝑖𝑥_𝑖(𝑘)−𝑥_𝑗(𝑘), (8) where 𝑎_𝑖𝑗>0 is a constant that expresses the effect level. According to Figure 2, we consider the following pairs of (𝑖,𝑗): ∈),(),(),(),(),(),(),(),(),(.𝑖,𝑗{(4,11,22,35,47,46,53,68,79,86,9)} (9) Other coupled terms 𝐴_𝑖𝑗 are given as 𝐴_𝑖𝑗=0. In addition, matrices 𝐴_𝑖𝑖,𝑖=1,2,…,9, can be obtained as follows: 𝐴₁₁=𝐴₁−𝐴₁₄,𝐴₂₂=𝐴₂−𝐴₂₁,𝐴₃₃=𝐴₃−𝐴₃₂,𝐴₄₄=𝐴₄−𝐴₄₅−𝐴₄₇,𝐴₅₅=𝐴₅−𝐴₅₆,𝐴₆₆=𝐴₆−𝐴₆₃,𝐴₇₇=𝐴₇−𝐴₇₈,𝐴₈₈=𝐴₈−𝐴₈₉,𝐴₉₉=𝐴₉−𝐴₉₆. (10) Thus, we can obtain the temperature model in multiple rooms as the form of (1).

Figure 1 

Illustration of multiple rooms.

Figure 2 

Directed graph expressing airflow.

3. Optimal Control Problem

For the large-scale system consisting of Σ_𝑖, 𝑖=1,2,…,𝑠, consider the optimal control problem. Define 𝑥(𝑘)∶=𝑥^𝑇₁(𝑘)𝑥^𝑇₂(𝑘)⋯𝑥^𝑇_𝑠(𝑘)^𝑇,𝑢_𝑖(𝑘)∶=𝑢^𝑐_𝑖(𝑘)^𝑇𝑢^𝑑_𝑖(𝑘)^𝑇𝑇,𝑢(𝑘)∶=𝑢^𝑇₁(𝑘)𝑢^𝑇₂(𝑘)⋯𝑢^𝑇_𝑠(𝑘)^𝑇. (11) Then, we consider the following problem.

Problem 1. Suppose that for the large-scale system consisting of subsystems Σ_𝑖,𝑖=1,2,…,𝑠, the current state 𝑥(𝑡)=𝑥_𝑡 is given, where 𝑡 is the current time. Then, find both continuous-valued and discrete-valued control inputs 𝑢(𝑘), 𝑘=𝑡,𝑡+1,…,𝑡+𝑁−1 minimizing the following cost function: 𝐽=^{𝑡+𝑁−1}_𝑘=𝑡𝑥(𝑘)−𝑥_𝑑^𝑇𝑄𝑥(𝑘)−𝑥_𝑑+𝑢(𝑘)^𝑇(𝑘+𝑥()𝑅𝑢)𝑡+𝑁−𝑥_𝑑^𝑇𝑄_𝑓𝑥(𝑡+𝑁)−𝑥_𝑑, (12) where 𝑄≥0, 𝑄_𝑓≥0, and 𝑅>0 are weighting matrices, and 𝑥_𝑑 is the desired state.

By assigning a binary variable to each element of 𝒰^𝑑_𝑖,𝑙, the subsystem Σ_𝑖 can be rewritten as a mixed logical dynamical (MLD) system [15], which is one of the standard models in hybrid systems. Therefore, Problem 1 can be rewritten as a mixed integer quadratic programming (MIQP) problem (see the appendix for further details). In the obtained MIQP problem, the dimension of continuous variables and that of binary variables are ∑^𝑠_𝑖=1𝑚^𝑐_𝑖𝑁 and ∑^𝑠_𝑖=1∑^𝑚𝑑𝑖_𝑙=1𝑝_𝑖,𝑙𝑁, respectively. On the other hand, the finite-time optimal control problem is frequently used in model predictive control (MPC). In MPC, the control input is generated as follows.

Procedure of MPC

Step 1. Set 𝑡=0 and give the current state 𝑥(𝑡)=𝑥_𝑡.

Step 2. Solve Problem 1.

Step 3. Apply only 𝑢(𝑡) to the plant.

Step 4. Set 𝑡∶=𝑡+1, measure 𝑥(𝑡), and return to Step 2.

In the above procedure, Problem 1, that is, the MIQP problem must be solved at each time. However, in large-scale systems, it is hard to solve Problem 1 within the practical computation time. Thus, it is necessary to consider a new method for solving Problem 1 under the situation that Problem 1 is used in MPC.

4. Proposed Method

In this section, first, the outline of the proposed method is explained. Next, the notion of virtual control inputs is proposed. Finally, using virtual control inputs, a solution method for Problem 1 is proposed.

4.1. Outline

For large-scale systems, it is in general difficult to solve Problem 1 within the practical computation time. Here, we focus on the fact that in MPC, only 𝑢(𝑡) is applied to the plant. From this fact, if the state at 𝑘=𝑡+1,𝑡+2,…,𝑡+𝑁 can be appropriately evaluated, then computation of the discrete-valued control input at 𝑘=𝑡+1,𝑡+2,…,𝑡+𝑁 is not required. In this paper, we consider relaxing the discrete-valued control input at 𝑘=𝑡,𝑡+1,𝑡+2,…,𝑡+𝑁 to a continuous variable, which is called here a virtual control input (see also Figure 3). In the proposed procedure of MPC, first, we find a continuous-valued control input and a virtual control input minimizing the cost function. Next, we derive only the discrete-valued control input at 𝑘=𝑡 from the virtual control input at 𝑘=𝑡.

Figure 3 

Illustration of state trajectory.

Furthermore, on implementation of the proposed procedure of MPC, we consider both centralized and decentralized controllers (see Figure 4). In the centralized controller, the continuous-valued control input and the virtual control input are computed. In each decentralized controller, the discrete-valued control input at 𝑘=𝑡 is computed using the virtual control input at 𝑘=𝑡. Thus, the computation load is decentralized.

Figure 4 

Control system considered in this paper.

Hereafter in this section, first, the notion of the virtual control input will be formally defined. Next, an approximate solution method of Problem 1 will be proposed. Finally, the proposed procedure of MPC will be shown.

4.2. Virtual Control Input

The matrix 𝐵^𝑑_𝑖 in subsystem (1) is rewritten as 𝐵^𝑑_𝑖=𝑆_𝑖𝐼_𝑟𝑖0𝑇00𝑇^1𝑖2_𝑖, (13) where 𝑟_𝑖∶=rank𝐵^𝑑_𝑖≤𝑚^𝑑_𝑖. Then, instead of subsystem (1), we consider the following subsystem: Σ^′_𝑖∶𝑥_𝑖(𝑘+1)=^𝑠_𝑗=1𝐴_𝑖𝑗𝑥_𝑗(𝑘)+𝐵^𝑐_𝑖𝑢^𝑐_𝑖(𝑘)+𝐵^𝑣_𝑖𝑢^𝑣_𝑖(𝑘), (14) where 𝐵^𝑣_𝑖=𝑆_𝑖𝐼_𝑟𝑖0. (15) The vector 𝑢^𝑣_𝑖(𝑘)∈𝒰^𝑣_𝑖⊆ℛ^𝑟𝑖 is the 𝑟_𝑖-dimensional continuous-valued control input. The set 𝒰^𝑣_𝑖 is given as 𝒰^𝑣_𝑖=[𝑢^𝑣_𝑖,𝑢^𝑣_𝑖]. The vector 𝑢^𝑣_𝑖 (𝑢^𝑣_𝑖) is given as the value of the vector 𝑇¹_𝑖𝑢^𝑑_𝑖(𝑘) such that a sum of each element of 𝑇¹_𝑖𝑢^𝑑_𝑖(𝑘) is minimized (maximized) and can be derived from the finite set 𝒰^𝑑_𝑖,𝑙.

Hereafter, 𝑢^𝑣_𝑖(𝑘) is called a virtual control input. We remark that the discrete-valued control input and the virtual control input are not one-to-one correspondence.

We show one example.

Example 2. Suppose that 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) in the subsystem (1) is given as 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘)=𝑢−1121−11^𝑑𝑖,1(𝑘)𝑢^𝑑_𝑖,2(𝑘)𝑢^𝑑_𝑖,3(𝑘), (16) where 𝑢^𝑑_𝑖,1(𝑘)∈{0,1,2}=∶𝒰^𝑑_𝑖,1,𝑢^𝑑_𝑖,2(𝑘)∈{0,2,4}=∶𝒰^𝑑_𝑖,2,𝑢^𝑑_𝑖,3(𝑘)∈{0,±3,±5,±7}=∶𝒰^𝑑_𝑖,3. (17) Noting here that rank 𝐵^𝑑_𝑖=2 holds, 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) in (16) can be rewritten as 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘)=−12𝑢11^𝑑𝑖,1(𝑘)−𝑢^𝑑_𝑖,2(𝑘)𝑢^𝑑_𝑖,3(𝑘). (18) Thus, we can obtain 𝐵^𝑣_𝑖𝑢^𝑣_𝑖(𝑘)=−12𝑢11^𝑣𝑖,1(𝑘)𝑢^𝑣_𝑖,2(𝑘). (19) In addition, the constraint for virtual control inputs 𝑢^𝑣_𝑖,1(𝑘) and 𝑢^𝑣_𝑖,2(𝑘) is given as −4−7_𝑢𝑣𝑖≤𝑢^𝑣_𝑖,1(𝑘)𝑢^𝑣_𝑖,2(𝑘)≤27_𝑢𝑣𝑖, (20) which can be computed from the finite sets 𝒰^𝑑_𝑖,1, 𝒰^𝑑_𝑖,2, and 𝒰^𝑑_𝑖,3.

Remark 3. Many actuators are included in each subsystem. In the example of air-conditioning systems, we can consider ceiling fans, local heaters, and so on. In many cases, these correspond to discrete-valued control inputs, and the number of these may be greater than the dimension of the state. Thus on the derivation of virtual control inputs, we reduce 𝑚^𝑑_𝑖 to 𝑟_𝑖. Another method is proposed in Section 5.

4.3. Solution Method for Optimal Control Problem

First, define 𝑢^′_𝑖(𝑘)∶=𝑢^𝑐_𝑖(𝑘)^𝑇𝑢^𝑣_𝑖(𝑘)^𝑇𝑇,𝑢^′(𝑘)∶=𝑢^′₁(𝑘)^𝑇𝑢^′₂(𝑘)^𝑇⋯𝑢^′_𝑠(𝑘)^𝑇𝑇. (21) Then, consider the following finite-time optimal control problem, instead of Problem 1.

Problem 2. Suppose that for the large-scale system consisting of subsystems Σ^′_𝑖,𝑖=1,2,…,𝑠, the current state 𝑥(𝑡)=𝑥_𝑡 is given, where 𝑡 is the current time. Then, find both continuous-valued and discrete-valued control inputs 𝑢^′(𝑘), 𝑘=𝑡,𝑡+1,…,𝑡+𝑁−1 minimizing the following cost function: 𝐽^′=^{𝑡+𝑁−1}_𝑘=𝑡𝑥(𝑘)−𝑥_𝑑^𝑇𝑄𝑥(𝑘)−𝑥_𝑑+𝑢^′(𝑘)^𝑇𝑅^′𝑢^′(𝑘+𝑥())𝑡+𝑁−𝑥_𝑑^𝑇𝑄_𝑓𝑥(𝑡+𝑁)−𝑥_𝑑, (22) where 𝑅^′ is the weighting matrix obtained from the weighting matrix 𝑅 in Problem 1 and the matrix 𝑇¹_𝑖 in (13).

Since the decision variable in Problem 2 is only a continuous-valued variable, Problem 2 can be rewritten as a quadratic programming (QP) problem, which can be solved by a suitable solver such as IBM ILOG CPLEX Optimizer [16]. The dimension of the decision variables is given by ∑^𝑠_𝑖=1(𝑚^𝑐_𝑖+𝑟_𝑖)𝑁. The procedure for rewriting Problem 2 as a QP problem is the same as that in the case of Problem 1. See the appendix for further details. In addition, from the definition of the virtual control input, we see that the minimum value of the cost function in Problem 2 gives the lower bound of the minimum value of the cost function in Problem 1.

By solving Problem 2, we can obtain the virtual control input. Next, we consider finding a discrete-valued control input from the obtained virtual control input. By 𝑢_𝑖^𝑣∗(𝑘), denote the optimal value of the virtual control input obtained by solving Problem 2. In addition, a binary variable 𝛿_{𝑖,𝑙,𝑞} is assigned to each candidate of each element in the discrete-valued control input. Then, focusing on only 𝑢_𝑖^𝑣∗(𝑡), consider the following problem for each subsystem Σ_𝑖.

Problem 3. Find binary variables 𝛿_{𝑖,𝑙,𝑞},𝑙=1,2,…,𝑚^𝑑_𝑖,𝑞=1,2,…,𝑝_{𝑖,𝑚𝑑𝑖} minimizing the following cost function: 𝐽_𝑑=𝑇¹_𝑖^𝑝𝑖,1_𝑞∑₌₁𝑢^𝑑_𝑖,1,𝑞𝛿_{𝑝𝑖,1,𝑞𝑖,2𝑞}∑₌₁𝑢^𝑑_𝑖,2,𝑞𝛿_{𝑖,2,𝑞𝑝𝑑𝑖𝑖,𝑚}⋮_𝑞∑₌₁𝑢^𝑑_{𝑖,𝑚𝑑𝑖,𝑞}𝛿_𝑖,𝑚^𝑑_𝑖,𝑞−𝑢_𝑖^𝑣∗(𝑡)₁, (23) under the equality constraint ∑^{𝑝𝑖,𝑙}_𝑞=1𝛿_{𝑖,𝑙,𝑞}=1, where ‖⋅‖₁ is the 1-norm of a vector.

By a simple calculation, Problem 3 can be rewritten as an integer linear programming (ILP) problem, which can be solved by a suitable solver such as IBM ILOG CPLEX Optimizer. We remark that Problem 3 is independent of each subsystem. That is, computing environment can be independently prepared for each subsystem. Furthermore, we remark that only the case of 𝑘=𝑡 is considered in Problem 3. This is because in MPC only the control input at 𝑘=𝑡 is applied to the plant. Thus, it is not necessary to find the discrete-valued control input at 𝑘=𝑡+1,𝑡+2,…,𝑡+𝑁−1. Therefore, Problem 3 can be solved fast. In addition, by deriving a look-up table off-line, solving the ILP problem online is not necessary. See the numerical example in Section 6 for further details.

The proposed solution method for Problem 1 can be regarded as the method that the MIQP problem corresponding to Problem 1 is divided into one QP problem (Problem 2) and 𝑠 ILP problems (Problem 3 for each subsystem). Problem 2 can be in general solved faster than Problem 1, and in Problem 3, we find only the discrete-valued control input at 𝑘=𝑡. From these observations, we see that the computation time for solving one QP problem and 𝑠 ILP problems is smaller than that for solving one MIQP problem.

In addition, the fact that Problem 1 is directly solved (in other words, Problem 1 is rewritten as an MIQP problem) corresponds to centralized control of large-scale systems. In the proposed method, solving Problem 2 (in other words, finding both continuous-valued control inputs and virtual control inputs) corresponds to centralized control. On the other hand, solving Problem 3 (in other words, finding discrete-valued control inputs based on optimal virtual control inputs) corresponds to decentralized control. Thus, the proposed method can be regarded as the method combining centralized control with decentralized control, and the computation load is decentralized and reduced.

4.4. Model Predictive Control Law

Finally, we summarize the procedure of MPC combining centralized control with decentralized control.

The Proposed Procedure of MPC

Step 1. Set 𝑡=0 and give the current state 𝑥(𝑡)=𝑥_𝑡.

Step 2. In the centralized controller, find both a continuous-valued control input and a virtual control input by solving Problem 2.

Step 3. Send the optimal values of both a continuous-valued control input and a virtual control input from the centralized controller to each decentralized controller.

Step 4. In each decentralized controller, solve Problem 3 and find a discrete-valued control input at 𝑡.

Step 5. Apply only the control input at 𝑡 to the plant.

Step 6. Set 𝑡∶=𝑡+1, measure 𝑥(𝑡), and return to Step 2.

5. Discussion on Quantization Errors

In the previous section, we do not consider quantization errors and focus on the lower bound of the cost function. In this section, we discuss a method for considering quantization errors.

First, the virtual control input is introduced. In this section, 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) in (1) is transformed into the virtual control input. By 𝑢^𝑞_𝑖(𝑘)∈ℛ^𝑛𝑖, denote the virtual control input. By 𝑢^𝑞_𝑖, denote the value of the vector 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) such that a sum of each element of 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) is minimized. By 𝑢^𝑞_𝑖, denote the value of the vector 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) such that a sum of each element of 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) is maximized. Then, instead of (1), we consider the following system: Σ_𝑖^′′∶𝑥_𝑖(𝑘+1)=^𝑠_𝑗=1𝐴_𝑖𝑗𝑥_𝑗(𝑘)+𝐵^𝑐_𝑖𝑢^𝑐_𝑖(𝑘)+𝑢^𝑞_𝑖(𝑘)+𝑎^𝑞_𝑖, (24) where 𝑢^𝑞_𝑖(𝑘)∈[𝑢^𝑞_𝑖,𝑢^𝑞_𝑖]⊆ℛ^𝑛𝑖. The vector 𝑎^𝑞_𝑖∈[𝑎^𝑞_𝑖,𝑎^𝑞_𝑖] is the uncertain affine term including quantization errors, and 𝑎^𝑞_𝑖,𝑎^𝑞_𝑖 are obtained from quantization errors in each element of 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘).

We show one example.

Example 4. Suppose that 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘) in subsystem (1) is given as 𝐵^𝑑_𝑖𝑢^𝑑_𝑖(𝑘)=𝑢110001^𝑑𝑖,1(𝑘)𝑢^𝑑_𝑖,2(𝑘)𝑢^𝑑_𝑖,3(𝑘), (25) where 𝑢^𝑑_𝑖,1(𝑘)∈{0,±1,…,±6}, 𝑢^𝑑_𝑖,2(𝑘)∈{0,±6}, and 𝑢^𝑑_𝑖,3(𝑘)∈{0,±1,…,±10}. Then, we can obtain the following virtual control inputs: 𝑢^𝑞_𝑖(𝑘)=𝑢^𝑑_𝑖,1(𝑘)+𝑢^𝑑_𝑖,2(𝑘)𝑢^𝑑_𝑖,3(𝑘). (26) The vectors 𝑢^𝑞_𝑖,𝑢^𝑞_𝑖 can be obtained as −12−10_𝑢𝑞𝑖≤𝑢^𝑞_𝑖(𝑘)≤1210_𝑢𝑞𝑖. (27) In addition, the following relations hold: 𝑢^𝑑_𝑖,1(𝑘)+𝑢^𝑑_𝑖,2(𝑘)∈{}𝑢−12,−11,…,12⊂𝒵,_𝑑𝑖,3(𝑘)∈{−10,−9,…,10}⊂𝒵. (28) From these relations, 𝑎^𝑞_𝑖,𝑎^𝑞_𝑖 can be obtained as 𝑎^𝑞_𝑖=,−0.5−0.5𝑎^𝑞_𝑖=0.50.5. (29) Thus, we can obtain virtual control inputs and the uncertain affine term.

Next, consider transforming the system (24) including uncertainty into a linear system including no uncertainty. For simplicity of discussion, we assume that, for all 𝑘≥0, two relations 𝑥_𝑖(𝑘)≥0 and 𝑢^𝑞_𝑖(𝑘)≥0 hold. This assumption can be satisfied by adding the affine term. Then, by using the result in [14], an overapproximation of the state can be expressed as 𝑥_𝑖(𝑘+1)∈𝑥_𝑖(𝑘+1),𝑥_𝑖(,𝑥𝑘+1)_𝑖(𝑘+1)=^𝑠_𝑗=1𝐴_𝑖𝑗𝑥^𝑐_𝑗(𝑘)−𝐴_𝑖𝑗𝑥^𝑟_𝑗(𝑘)+𝐵^𝑐_𝑖𝑢^𝑐_𝑖(𝑘)+𝑢^𝑞_𝑖(𝑘)+𝑎^𝑞_𝑖+𝑏,𝑥_𝑖(𝑘+1)=^𝑠_𝑗=1𝐴_𝑖𝑗𝑥^𝑐_𝑗(𝑘)+𝐴_𝑖𝑗𝑥^𝑟_𝑗(𝑘)+𝐵^𝑐_𝑖𝑢^𝑐_𝑖(𝑘)+𝑢^𝑞_𝑖(𝑘)+𝑎^𝑞_𝑖+𝑏, (30)(31) where 𝑏 is the affine term for satisfying the above assumption, and 𝑥^𝑐_𝑖(𝑘)𝑥^𝑟_𝑖(𝑘)=12𝐼_𝑛𝐼_𝑛−𝐼_𝑛𝐼_𝑛𝑥_𝑖(𝑘)𝑥_𝑖(𝑘). (32) Equations (30) and (31) are a discrete-time linear system, and by using (31), the worst case of the cost function (12) in Problem 1 can be minimized approximately. This problem is reduced to a quadratic programming (QP) problem, not a min-max optimization problem. On the other hand, by using both (30) and (31), the middle point of the cost function (12) can be evaluated. Also in this case, the optimal control problem is reduced to a QP problem.

Remark 5. In this section, the dimension of the virtual control input is 𝑛_𝑖. On the other hand, in Section 4.2, the dimension of the virtual control input is 𝑟_𝑖=rank𝐵^𝑑_𝑖≤𝑛_𝑖. Also in the discussion on quantization errors, the dimension of the virtual control input can be reduced to 𝑟_𝑖. However, 𝑢^𝑞_𝑖,𝑢^𝑞_𝑖 and 𝑎^𝑞_𝑖,𝑎^𝑞_𝑖 must be determined depending on the choice of 𝑇¹_𝑖,𝑇²_𝑖 in (13). In order to avoid this complexity, we have considered in this section the virtual control input with the dimension of 𝑛_𝑖.

6. Numerical Example

As a numerical example, consider the large-scale system consisting of 9 subsystems given by (1). Couplings of subsystems are given by the directed graph in Figure 5. The matrices 𝐴_𝑖𝑖, 𝑖=1,2,…,9, are given by 𝐴_𝑖𝑖=𝑎_𝑖0.601, (33) where 𝑎₁=1.3, 𝑎₂=1.2, 𝑎₃=1.1, 𝑎₄=1.0, 𝑎₅=0.9, 𝑎₆=0.8, 𝑎₇=0.7, 𝑎₈=0.6, and 𝑎₉=0.5. The matrices 𝐴_𝑖𝑗, 𝑖≠𝑗 are given as follows. 𝐴₁₂, 𝐴₂₁, 𝐴₂₃, 𝐴₃₂, 𝐴₁₄, 𝐴₄₁, 𝐴₂₅, 𝐴₅₂, 𝐴₃₆, 𝐴₆₃, 𝐴₄₅, 𝐴₅₄, 𝐴₅₆, 𝐴₆₅, 𝐴₄₇, 𝐴₇₄, 𝐴₅₈, 𝐴₈₅, 𝐴₆₉, 𝐴₉₆, 𝐴₇₈, 𝐴₈₇, 𝐴₈₉, 𝐴₉₈ are given as 0.4𝐼₂, respectively. Other coupling terms 𝐴_𝑖𝑗, 𝑖≠𝑗 are given as a zero matrix. 𝐵^𝑐_𝑖 is given by 𝐵^𝑐_𝑖=01. (34) In addition, we impose the following constraint: −10≤𝑢^𝑐_𝑖(𝑘)≤+10. (35) The matrix 𝐵^𝑑_𝑖 and the discrete-valued control input 𝑢^𝑑_𝑖(𝑘) are given by those in Example 2.

Figure 5 

Directed graph expressing couplings.

For this system, we consider solving Problem 2, instead of Problem 1. The prediction horizon, the weighting matrices, the desired state, the initial time, and the initial state are given as 𝑁=5, 𝑄=10𝐼, 𝑅^′=𝐼, 𝑥_𝑑=0, 𝑡=0, and 𝑥_𝑖(0)=[108]^𝑇, respectively.

We show the computation result. Figures 6 and 7 show state trajectories of 9 subsystems. Note here that in these figures, the virtual control input is used. From this fact, we see that trajectories in these figures imply those corresponding to the lower bound of the optimal value of a given cost function.

Figure 6 

Trajectories of the first element of the state.

Figure 7 

Trajectories of the second element of the state.

Next, we discuss how to derive the discrete-valued control input from the virtual control input. In this example, the optimal virtual control input at 𝑘=0 can be obtained by 𝑢_𝑖^𝑣∗(0)=2−7,𝑖=1,2,…,9. (36) Then focusing on (18), we can obtain the discrete-valued control input as follows: 𝑢^𝑑_𝑖(0)=𝑢^𝑑_𝑖,1(0)𝑢^𝑑_𝑖,2(0)𝑢^𝑑_𝑖,3(0)∈20,42−7−7. (37) In this case, 𝑢^𝑑_𝑖,3(0) can be uniquely determined. 𝑢^𝑑_𝑖,1(0) and 𝑢^𝑑_𝑖,2(0) are not uniquely determined. If it is desirable to set a small value, then we can obtain 𝑢^𝑑_𝑖,1(0)=2 and 𝑢^𝑑_𝑖,2(0)=0. In practice, the look-up table can be generated off-line. Then, the input of this table is the rounded virtual control input, and the output is the discrete-valued control input. In online computation, the discrete-valued control input can be determined from the look-up table and the rounded virtual control input. In this method, it is not necessary to solve Problem 3.

Finally, we explain the computation time for solving Problem 2. In this example, the computation time of the QP problem was 0.05 [sec], where we used IBM ILOG CPLEX Optimizer 11.0 as an MIQP solver on the computer with the Intel Core 2 Duo 3.0 GHz processor and the 4 GB memory. Thus, Problem 2 can be solved fast.

7. Conclusion

In this paper, we have proposed computational techniques for model predictive control of large-scale systems with both continuous-valued and discrete-valued control inputs. By introducing the notion of virtual control inputs, the computation load is reduced and decentralized. The effectiveness of the proposed method has been shown by a numerical example. The proposed method is useful for solving the control problem of large-scale system with both continuous-valued and discrete-valued control inputs.

One of the future works is to apply our approach to several practical systems such as air-conditioning systems. On the other hand, a large-scale system studied in this paper is a class of hybrid systems. It is also important to extend the proposed method to general hybrid systems such as piecewise affine systems and MLD systems.

Appendix

Reduction of Problem 1 to an MIQP Problem

First, by using a binary variable, the element 𝑢^𝑑_𝑖,𝑙(𝑘) of 𝑢^𝑑_𝑖(𝑘) in subsystem (1) can be expressed as 𝑢^𝑑_𝑖,𝑙(𝑘)=𝑢^𝑑_𝑖,𝑙,1⋯𝑢^𝑑_{𝑖,𝑙,𝑝𝑖,𝑙}𝛿_𝑖,𝑙,1(𝑘)⋮𝛿_{𝑖,𝑙,𝑝𝑖,𝑙}(𝑘)=𝐶^𝑑_𝑖,𝑙𝛿_𝑖,𝑙(𝑘), (A.1) where we imposed the following equality constraint: 𝛿_𝑖,𝑙,1(𝑘)+⋯+𝛿_{𝑖,𝑙,𝑝𝑖,𝑙}(𝑘)=1. (A.2) By using 𝐶^𝑑_𝑖,𝑙 and 𝛿_𝑖,𝑙(𝑘), 𝑢^𝑑_𝑖(𝑘) can be expressed as 𝑢^𝑑_𝑖(𝑘)=𝐶^𝑑_𝑖,10⋱0𝐶^𝑑_{𝑖,𝑚𝑑𝑖}𝛿_𝑖,1(𝑘)⋮𝛿_𝑖,𝑚^𝑑_𝑖(𝑘)=𝐶^𝑑_𝑖𝛿_𝑖(𝑘). (A.3) Next, by combing all subsystems, we can obtain 𝑥(𝑘+1)=𝐴𝑥(𝑘)+𝐵𝑢(𝑘), (A.4) where 𝐴𝐴=¹¹⋯𝐴_1𝑠𝐴⋮⋱⋮_𝑠1⋯𝐴_𝑠𝑠,𝐵𝐵=^𝑐1𝐵^𝑑₁𝐶^𝑑₁0⋱0𝐵^𝑐_𝑠𝐵^𝑑_𝑠𝐶^𝑑_𝑠. (A.5) See Section 3 for the definitions of 𝑥(𝑘)and𝑢(𝑘). In addition, the constraint for the continuous-valued control input and the equality constraint of (A.2) can be expressed as 𝐶𝑥(𝑘)+𝐷𝑢(𝑘)≤𝐸. (A.6) Since 𝑢(𝑘) consist of both continuous variables and binary variables, the pair of (A.4) and (A.6) is a class of mixed logical dynamical systems [15]. From (A.4) and (A.6), we have 𝑥=𝐴𝑥₀+𝐵𝑢,𝐶𝑥+𝐷𝑢≤𝐸, (A.7)(A.8) where 𝑥∶=𝑥^𝑇(𝑡)𝑥^𝑇(𝑡+1)⋯𝑥^𝑇(𝑡+𝑁)^𝑇,𝑢∶=𝑢^𝑇(𝑡)𝑢^𝑇(𝑡+1)⋯𝑢^𝑇(𝑡+𝑁−1)^𝑇,𝐼𝐴𝐴=𝐴²⋮𝐴^𝑁,𝐴𝐵=00⋯0𝐵0⋯0𝐴𝐵⋱⋱⋮⋮⋱⋱0_𝑁−1,𝐵⋯𝐴𝐵𝐵𝐶=,𝐶00⋱⋮0𝐶0𝐷=⋱,𝐷00𝐷𝐸⋮𝐸𝐸=. (A.9) By substituting (A.7) into (A.8), we have 𝐷−𝐶𝐵𝑢≤𝐸−𝐶𝐴𝑥₀. (A.10) Finally, consider the cost function (12). For simplicity of discussion, we consider the case of 𝑥_𝑑=0. Then the cost function (12) can be rewritten as 𝐽=𝑥^𝑇𝑄𝑥+𝑢^𝑇𝑅𝑢, (A.11) where ⋱𝑄𝑄=𝑄00𝑄_𝑓,𝑅=⋱𝑅00𝑅. (A.12) Furthermore, by substituting (A.7) into this cost function, we can obtain 𝐽=𝑢^𝑇𝑅+𝐵^𝑇𝑄𝐵𝑢+2𝑥^𝑇₀𝐴^𝑇𝑄𝐵𝑢+𝑥^𝑇₀𝐴^𝑇𝑄𝐴𝑥₀. (A.13) From the above, Problem 1 can be rewritten as the following MIQP problem: find a control input sequence 𝑢 minimizing the cost function (A.13) subject to the constraint (A.10). The obtained MIQP problem can be solved by a suitable solver such as IBM ILOG CPLEX Optimizer. However, for large-scale systems, it is difficult to solve the MIQP problem with practical computation time.

Acknowledgment

This work was partially supported by Grant-in-Aid for Young Scientists (B) 23760387.