Model predictive control-based scheduler for repetitive discrete event systems with capacity constraints

. A model predictive control-based scheduler for a class of discrete event systems is designed and developed. We focus on repetitive, multiple-input, multiple-output, and directed acyclic graph structured systems on which capacity constraints can be imposed. The target system’s behaviour is described by linear equations in max-plus algebra, referred to as state-space representation. Assuming that the system’s performance can be improved by paying additional cost, we adjust the system parameters and determine control inputs for which the reference output signals can be observed. The main contribution of this research is twofold, 1: For systems with capacity constraints, we derived an output prediction equation as functions of adjustable variables in a recursive form, 2: Regarding the construct for the system’s representation, we improved the structure to accomplish general operations which are essential for adjusting the system parameters. The result of numerical simulation in a later section demonstrates the effectiveness of the developed controller.


Introduction
The objective of this research is to design and implement a model predictive control (MPC) [1,2] based scheduler for a class of discrete event systems (DESs) [3,4].A discrete event system is a dynamical system whose state evolves at discrete time instants.
The dynamics of such systems can be formulated by a pair of linear equations in maxplus algebra [8,9], and a system whose behaviour is described by these equations is referred to as a max-plus linear system.Typical examples include manufacturing and transportation systems [10,11].Since the form of the pair is similar to the state-space representation in modern control theory, several concepts in control theory have been applied to max-plus linear systems [12][13][14][15].Amongst the related concepts, the frameworks for controlling manufacturing systems based on MPC provide cost effective schedulers, which can be roughly classified into two types.
One type deals with systems whose feature quantities of the system parameters are given [6], [16,18].The control objective is to determine an optimal input time that minimises an objective function for given reference output times; the function includes an output error and a rate of change of input times.References [6], [16] consider a deterministic system, while [17] extends this to a stochastic case.Reference [18] takes into account external perturbation and considers stabilizing a system by feedback control.A typical feature of these works is that the system parameters are either constant or vary stochastically around predetermined values.
The other type of framework handles systems in which the system parameters can be adjusted as well as the input times.For example, the framework in reference [7] first adjusts the system parameters to minimise both an operational cost and tardiness penalty, and then determines an input time which minimises a waiting cost.If we apply the framework to manufacturing systems, the controller aims to reduce the manufacturing cost, avoid delays for the due dates, and supply materials at as late as possible.
Reference [19] proposed an efficient method to handle the system parameters in the output prediction equation as adjustable variables by introducing the concept of a 'cell'.The control framework is reduced to solving two optimisation problems in sequence.However, the application scope is currently restricted; the framework assumes that only a single job can be processed at one time in a single facility, and that the number of in-process jobs, which can exist between two adjacent facilities, is unlimited.
However, capacity constraints are occasionally incurred on several locations in most actual systems.Thus, a framework to control such systems needs be developed.In view of this, we develop an MPC based scheduler in which maximum capacity constraints in a single or between two arbitrary facilities can be incurred.The main contribution of this research is twofold.
The first contribution is regarding the output prediction equation.Starting from the state-space representation for systems with capacity constraints, we derived a prediction equation as functions of adjustable variables.Since it is hard to represent the resulting matrices equation in an explicit form, each block of the matrices is expressed in a recursive form.
The second contribution is an improvement of the structure of the 'cell'.The previous structure is designed only for two special operations.However, in handling systems with capacity constraints, a framework for general addition and multiplication is required to construct an optimisation problem for adjusting the system parameters.We thus improve the structure and create procedures for these operations.

Algebraic system
We first define a set , where R represents the real set.Then, for x , y  max R , we define three operators and two unit elements: . In a conventional formalism, the special rule , and , we define

, and
. For the unit matrices, ε is a matrix, all elements of which are  , while E is a matrix whose diagonal elements are e and off-diagonal elements are  .The priority of operators  and \ is higher than that of  and  , and  is often omitted when no confusion is likely to arise.If and b N  , we define and write the b -th power of X by

State-space representation
Let the token number, number of nodes, and number of external inputs and outputs be k , n , r , and m , respectively.We denote the processing, process start, process completion, input, and output times for the k -th token by Moreover, we introduce the following representation matrices: H e : maximum number of tokens between node i and its downstream node j is h .
, and are referred to as the weight, adjacency, capacity, input, and output matrices, respectively.
Then, the earliest process start, process completion, and output times of the k -th token are represented by the following equations [20]: where 1) and ( 2) are known as the state and output equations, respectively.
Q is the maximum capacity of the system.Operator * is called the Kleene star.For a DAG structured adjacency matrix

Design of an MPC controller
The design of an MPC controller generally consists of the following three steps.1) Derive prediction and output equations 2) Determine control inputs 3) Apply the receding horizon method In the first step, we apply the state equation Eq. ( 1) iteratively and derive prediction equations for the tokens, where N is known as the prediction horizon.We then apply the output equation Eq. ( 2) to the prediction equations, for which the output times of the corresponding tokens are predicted.
In the second step, we set reference signals for the output equations, which correspond to give target output times.We then construct an optimisation problem to determine a control law.To reduce the computation load, a technique to decrease the number of variables is often used.Given a certain constant , we assume that the adjustable variables for the


-th tokens remain unchanged, where c N is known as control horizon.The variables include the control input ) (k u and system parameter ) (k d .In respect to the optimal parameters for the -th tokens, we apply only the first one, and the succeeding laws are determined after the token number is incremented and the parameters recalculated.This method is known as the receding horizon method, by which a feedback control is achieved.

Proposed Framework
In developing an MPC framework for systems with capacity constraints, first, we derive the prediction and output equations.Then, the system parameters are adjusted by solving an optimisation problem, and subsequently the optimal input time is determined.

Prediction equation
The prediction equation can be derived using Eq.(1).Letting the prediction horizon be N , we formally express the prediction equation for the token as follows: For 0  l , by comparing the coefficients of 1) and ( 4), we have If 1  l , we advance the token number k in Eq. ( 1) for l steps, by which is obtained.On the other hand, by replacing l with m l  in Eq. ( 4), we obtain If Q l  , we substitute this for the first term on the right-hand side in Eq. ( 6).This yields The comparison of the coefficient for 4) and ( 8) yields On the other hand, if Q l  , Eq. ( 6) can be expanded as The second term on the right-hand side can be transformed into Using Eqs. ( 7) and ( 11), Eq. ( 10) is represented as follows: Let us compare the coefficients of ) ( h k  x in Eqs. ( 4) and ( 12 , we obtain the following relationship On the other hand, if , the third term on the right-hand side in Eq. ( 12) can be omitted and we obtain In view of Eqs. ( 9) and ( 14), h l , Φ can be expressed in the following unified manner 4), (8), and ( 12).If l p  , only the last terms in these equations are effective and we obtain 4) and ( 8) yields 4) and (12) Equations ( 16) and ( 17) can be then unified to for l p  .Furthermore, we derive an output prediction equation for the 2) and ( 4), the simultaneous output prediction equations for these tokens can be expressed as where In this manner, the coefficient matrices Φ and Δ in the output prediction equation ( 19) are now derived.These matrices are functions of the system parameters N l handled as adjustable variables, each block of which is expressed in a recursive form.

Control objective
The set of output prediction equations expressed in Eq. ( 19) is used to construct a control law for which the given reference output times are observed.
Let the reference output times for the , respectively.The well-known control laws are to adjust where We assume that the operational and waiting costs can be reduced if we can enlarge , respectively.On the other hand, a tardiness penalty is incurred if an actual output time exceeds the reference time.A control law is then designed by taking into account these three points.Let l μ be the tardiness vector for the One natural idea for control objective is, recalling Eq. ( 21), to enlarge , and decrease where

R
can be interpreted as the set of artificial variables.
However, if we formulate this objective in a straightforward manner, the number of variables and constraints increases very sharply as the prediction horizon increases, and the computation time inflates accordingly.Hence, we use a lightweight method used in [21] which is understood as a greedy method; we determine in two separate steps.

Adjustment of the system parameters
An optimisation problem is constructed and solved to adjust the system parameters . First, we ignore the effect of the input times ) (k U in Eq. ( 22), and formulate a control objective as where Recalling Eq. ( 20), the j -th row of the l -th row block in Eq. ( 23) can be rewritten as . This can be further reduced to the following simultaneous inequalities.
for all h , i , j , and l ( . For the other constraints, we assume that the lower and upper bounds on the system parameters are as follows: With respect to the objective function, we consider reducing both the operational cost and tardiness, and formulate the objective function as Since this problem has the form of a linear programming problem, we can solve the problem by using an efficient algorithm such as the interior point method [22].
In addition to these, we consider introducing the concept of control horizon to reduce the number of variables.Letting the control horizon be c N , we assume Applying this concept, the number of system parameters can be decreased from n N  to n N c  .According to [21], we can obtain sufficient performance even if we set 1  c N .

Determination of the input times
In the second step of the control law, the input times ) (k U are determined.If the control inputs are small enough not to affect the output times, the earliest output times are predicted as If the values of the state vector are sufficiently large, ) (k Y may exceed the reference output time ) (k R .For such cases, we aim to adjust the input times to complete the k -th token at ) (k Y and not at ) (k R .This policy can be expressed as We should note here that the term ) (k X Φ appears on the right-hand side compared to Eq. (21).In Eq. (29), only ) (k U is a variable vector, while the others are constant vectors and matrices.
To reduce the waiting cost, the input times should be as large as possible in observance of Eq. ( 29).Using a lightweight formula known as the residuation theory [5], the optimal solution for ) (k U is calculated as follows: The optimal input times ) are then obtained using this result.In the last step, we consider applying the concept of receding horizon method to achieve a feedback control.For the obtained optimal input times, we give an actual control input only for the first step.The inputs for the subsequent steps are recalculated after k is incremented and the values of the latest state vector are obtained.We adopt and give only the first block of Eq. (30) for the optimal input, the specific form of which is expressed as where Using the optimal solutions given in Eqs. ( 28) and (31), we can operate the system at a lower operational cost observing the reference output times.

Improvement of the cell's structure
To solve the optimisation problem constructed in section 3.3, a method for deriving and representing h l, Φ in Eq. ( 25) as a function of the system parameters is essential.This matrix is calculated iteratively using Eqs.( 5), ( 13), (15), and ( 18), and k A includes the Kleene star of the weighted adjacency matrices as shown in Eq. ( 3).
Reference [19] proposed a method for computing the transition matrix in the state-space representation as a function of the system parameters by introducing the concept of a 'cell'.However, the method is designed only for calculating the Kleene star and multiplying a diagonal matrix.
On the other hand, in handling systems with capacity constraints, a framework for general addition and multiplication in max-plus algebra is required.Thus, we improve the structure of a cell.
Figure 1 depicts the improved structure and an example.Each resulting matrix is described by collection of cells and each element of which is described by one cell.For a given matrix X , Src and Dst store the column and row numbers of matrix X , respectively.If there is no cell that has (Src, Dst)= ) , ( i j , this indicates . Arrays Coef() and Pows() are used for representing the coefficient and power of a function of the system parameters . Arrays Mins() and Maxs() are used for recording the minimum and maximum weighted sums for given lower and upper bounds of the system parameters We now consider the case for X

d
, and 2   c N N .Then, the cell in Figure 1 indicates [d 6.Using this structure, the resulting matrices in the state-space representation and output equations can be represented as functions of the system parameters.Since the computations for the Kleene star of a matrix and the multiplication of a diagonal matrix have already been developed in [19], we newly develop the remaining operations required to derive the output prediction equation.

Addition
Two given cells with equal Src and Dst can simply be added by concatenating the Coef(), Pows(), Mins(), and Maxs() of the two cells.However, a naïve concatenation of these may yield redundant terms.Thus, a procedure for removing the redundant terms, which is considered as a branch and bound [23], is needed.This procedure is repeated for all rows and columns of the given matrices.

Multiplication
To compute Z X  for matrices and , the resulting cell is equal to cell z and (Src, Dst) is changed to ) , ( i j , and vice versa.The addition of the resulting cell to the previous result is repeated for all l ) 1 ( n l   .

Branch and bound
After the above addition or multiplication above, the resulting cell may have redundant terms that have no effect on the result.Such terms should be removed from the cell to reduce the computation time.For a target cell x, we can remove the terms that cannot be the maximum value.In the first step, we remove terms s In the case of Figure 1, holds and the second term 2  s is removed.Next, we evaluate the order relations for all pairs of terms i and j  26), the significant number of resulting terms would be reduced.In view of this, this procedure can be seen as a kind of branch and bound.

Difference from the previous structure
The difference between the previous structure in [19] and the proposed one in this paper is that element Dst, and arrays Mins() and Coef() are appended.By appending Dst, the structure of a cell becomes self-descriptive and general operation for multiplication is achieved.Array Mins() is used in the branch and bound procedure, and plays a role in reducing the complexity of the resulting matrix.The original objective of appending array Coef() is to improve the descriptive ability, and the direct benefit here is that a lag time between two successive processes can be taken into account.

Numerical Simulation
We present an application example of a manufacturing system.Figure 2  external output.The values in parentheses (*) represent the system parameters, which correspond to the processing times, while the values in brackets [*] represent the maximum number of in-process jobs that can exist in a single facility or between two facilities.For instance, facility 2 can process a maximum of three jobs, and the maximum number of inprocess jobs that can exist between facilities 2 and 3 is five.Since the ratio of the maximum capacity [*] to the system parameter (*) gives the effective throughput (jobs per unit time) for each facility, facility 4 provides the highest throughput.We assume that the transition times between two adjacent facilities, starting upstream, are 0.2, 0.3, and 0.4, respectively.By incurring an additional operational cost, the processing times can be reduced, the minimum values of which are 0.5 less than the original ones.Then, we have the following representation matrices and vectors:

H
are generated by executing the following computation for all h For the weights in the objective function Eq. ( 27), we assume:  for varying prediction horizons N =1, 2, 5, 10, and 15 with control horizon N =1.At the beginning stage, the parameter is set to the normal value (=5.0) for all cases.Then, for N =1, it decreases sharply to the minimum value 4.5 at  k 9 and 23.This phenomenon occurs 2 or 3 steps after the interval of the reference output time has been decreased.This lag is due to the capacity of facility 3.For N =2, the parameter is decreased at an earlier stage than for N =1.For N =5 and 10, the parameter decreases and recovers more smoothly at even an earlier stage, and does not decrease to the minimum limit.Furthermore, for N =15, the parameter does not decrease and remains constant all through the jobs.
Figure 4 depicts the output error for N =1, 2, 5, 10, and 15 with c N =1.We define the error for job , which can be interpreted as tardiness for the reference output time.For N =1 and 2, tardiness occurs soon after the interval of the reference output time has been decreased.By contrast, for N =5, the output time is moved up before  k 5 and 20 to avoid tardiness, but tardiness eventually occurs at  k 10 and 25.For N =10, the output times are moved up at 9 1   k and 24 12   k , and tardiness is completely avoided.Tardiness is also avoided for N =15 and the output times are further moved up.In the light of Figures 3 and 4, if we set the prediction horizon N to a large value, we can avoid tardiness for the reference output time with only a small additional operational cost.
Table 1 shows the number of constraints for the linear programming problem in Eq. ( 28) and the performance of the controller.We define the performance index T as where the first and second terms represent the sum of the operational cost and penalty for tardiness, respectively.Clearly, the performance is improved as the prediction horizon N increases.On the other hand, since the number of constraints increases, so too does the computation time to obtain the internal representation of the output prediction equations and to solve the linear programming problem.

Conclusion
We have developed an MPC controller for a class of MIMO-FIFO discrete event systems with adjustable system parameters, in which capacity constraints exist in a single facility or between two arbitrary facilities.We first derived an output prediction equation as functions of the system parameters, each block of the resulting matrices of which is obtained in a recursive form.
We then improved the structure of a cell for system's internal representation, and created procedures for general addition and multiplication in max-plus algebra.Through a numerical simulation, we have demonstrated that the proposed MPC controller achieves good performance if the prediction horizon is set to a large value.
for operational cost and tardiness, respectively.Consequently, the optimisation problem for ) ( l k  d and l μ is summarised as min  J subject to Eqs. (24)-(26).
arrays' size is equal to the number of irreducible terms.Let S and s be the number of irreducible terms and the term number in the cell, respectively.Then, the type and size of the arrays are Coef()

Figure 1 .
Figure 1.The structure and an example of a cell.
For two given cells x and z that have (Src, Dst).Then, recalling the property of the distributive law, arrays Coef(), Pows(), Maxs(), and Mins() are multiplied (added in conventional algebra) for all pairs of of term i is constantly equal to or smaller than term j , and thus term i can be removed from the current cell.If tight bounds were set for ) depicts a simple tandem structured manufacturing system with four facilities, one external input, and one

Figure 3
Figure 3 illustrates the profile of the system parameter ) ( 2 k d for varying prediction horizons

Figure 2 .Figure 3 .
Figure 2. A simple system with four facilities

Table 1 .
The number of constraint inequalities.
[21] Goto, H., A lightweight model predictive controller for repetitive discrete event systems, Asian Journal of Control, 2013 (in press) [22] Nocedal, J. and Wright, S., Numerical Optimization Second Edition.Springer, New York (1999).[23] Hillier, F. S. and Lieberman, G. J., Introduction to Operations Research.McGraw Hill, New York (2004).Hiroyuki Goto is a professor at the department of industrial & system engineering, Hosei University, Japan.He received his B.S. and M.S. degrees from The University of Tokyo in 1995 and 1997, respectively.He received his D.E. degree from Tokyo Metropolitan Institute of Technology in 2004.His research interests include controller design for discrete event systems and high-performance computing.