Design of a Power System Supervisory Control with Linear Optimization for Electrical Load Management in an Aircraft On-Board DC Microgrid

In modern aircraft, energy supply management has become a critical matter, since many aboard electrical loads have to be supplied, especially those related to flight safety. However, at the same time, the size and weight of electrical generators must be limited because of their on-board installation. In this paper, the Mixed Integrated Linear Programming (MILP) methodology has been used to formulate the Supervisor definition of the direct current (DC) microgrid (MG) on-board system with an extension for the programmable loads. Due to the problem of dimension increase, two methods have been presented and tested to perform optimal energy management (EM) aboard an aircraft: the Branch and Bound (B&B) and the Linear Regression Approximation (LRA). Finally, numerical simulations and results have been provided to validate the proposed optimization methodologies, according to the dimensions and the complexity of the problem.


Introduction
Over the last few years, the electrical power distribution aboard an aircraft has become a critical issue because of the rising number of electrical devices installed on-board, which leads to an increase in supply demand. On the other hand, the size of the power generators has to be limited in order to reduce both their weight and occupied volume, at the same time preventing any shortage of energy in order to maintain the correct operation mode of all the on-board supplied loads.
The power limits impose the rigorous control of the power itself; for this reason, the on-board power grid commonly integrates protection circuits, remote controllable circuit breakers and power converters with defined regulators and logics, in order to improve the power quality and, at the same time, to avoid generator overload. The logics have the task of disconnecting the low-priority loads, in the case of generator overload, and of reconnecting them when the overload is over, or, otherwise, controlling their power de-rating.
The family of logics and controls implemented on-board for optimal power flow, including storage and sources, is generally called Energy Management (EM) and it runs in a single or distributed board computers, commonly called EM Supervisors. Additional terms are also introduced to better explain the EM domain: for example, the logics that have loads as their domain are called Electrical Load Management (ELM), while, if the logic works at the source level, we simply refer to the terms Source Management. Finally, if the grid has low storage, we commonly refer to Power Management (PM). In this paper, we generally refer to the EM term for the logics and control methods, since it is the most general; reference is made to Supervisor to generally indicate the on-board computer or the implementation of the algorithm itself, where the EM logics are performed [1][2][3].
ing and reconnecting loads according to their importance, in addition to power regulation, has many advantages and disadvantages, illustrated in various papers [1,20].
EPC is an Electrical Power Center for aeronautical applications, which deals with the power regulation of the connected loads; it is equipped with solid-state power controllers (SSPC), which communicate via a Controller Area Network (CAN or CAN Bus) type 2.0B with a supervisor, where the energy management logic is executed [3]. The control is able to disconnect or reduce the power to non-essential loads automatically, to prevent the main generator from being overloaded [21].
The EPC requires two continuous-voltage inputs: 270 V DC for power loads supply and 28 V DC for feeding the digital parts (micro-controllers, communication devices, drivers). The details of the required power for each load are shown in Table 1, comprising also the two electromechanical actuators (EMA), which are the Landing Gear System (LGS) and Flight Control System (FCS), besides the other connected passive loads (R1-R5): The EPC configuration of the flight test stand connections is shown in Figure 1.   (ECS) supplies an alternating current of 115 V AC with variable frequency between 400 and 800 Hz; the transformerrectifier unit (TRU1, TRU2) levels up to a 270 V DC output voltage to feed the main power  bus of the EPC; a DC/DC converter provides the 28 V DC voltage from the 270 V DC voltage bus previously obtained, to supply both EPCs and the other electrical devices aboard the aircraft.

Description of a Microgrid
Current EPC systems must have a power management strategy to control the flow of power in a microgrid system. A typical microgrid is composed of the AC voltage distribution grid interfaced with a VSI (Voltage Source Inverter) [22], which acts as a grid formation converter, an energy storage system consisting of a battery bank, and a bidirectional DC-DC converter [23].
In any DC microgrid system (MG), it is important to use control strategies that include both stability analysis and cell stabilization techniques [24]. Control is usually classified into two levels: local and coordinated. The local one is based only on local measurements and only through a communication line can it be made available among the various microgrids in order to achieve coordinated control. From the method of communication, three strategies can be obtained: basic coordinates (or decentralized control), centralized and distributed [25,26]: decentralized control [27] can be considered an extension of local control and it is based only on local measurements, while centralized and distributed control strategies are based on digital communication technologies [28,29].
The correct distribution of resources, while maintaining system efficiency and reliability, is one of the main priorities of a DC microgrid system. Indeed, for each type of DC-DC converter, there are many advantages and disadvantages to be taken into consideration with regard to the hierarchical control strategies: primary and secondary [30]. A bidirectional regulation system can consist of a transformer-less H-bridge inverter used to regulate the power flow and to limit the fault current flowing in the distribution grid: the inverter control strategy checks the power flow at the common coupling point (PCC), while the operation of the H-bridge allows us to limit the short-circuit current (or, at least, to regulate it) during a disturbance in the network, fostering its reinstatement after the fault has been eliminated [31]. Furthermore, accurate transient modeling of faulty DC microgrids is the basis of fault detection location and isolation [32].
The power management technique applied to the concept of virtual inertia together [33] with a power management function of the generators allows us to avoid high power peaks [34] and improves the availability of energy by preventing and without degrading the DC bus voltage.

Electrical Load Management Definition by MILP
In this section, the problem is formulated as Mixed integer Linear Programming (MILP); the core of the control scheme is reported in Figure 2. For our discussion, we consider the plant status as fully observable, thus allowing the calculation, for each time-step, of the available power P AV . The MILP based power management control has additional inputs, which are the load priorities, their power size and the vector ν, having the elements set to 1, in correspondence with the loads that can be reduced in power, or to 0 for ON/OFF loads.
The elements of the output vector x can assume 0 value for the OFF state, 1 value for the ON state and a number between 0 and 1 for those loads allowing the power reduction.

Assumptions
In this paper, we consider only the optimal steady-state problem solution; this implies that the transient is not studied and all the variables are linearly related to each other. The Thevenin's Theorem is adopted to model the generator. Moreover, as a simplifying hypothesis, the grid is fully observable, i.e., the power of the generator and that of the loads are known.
As an additional assumption, we consider active loads, with power reference as input. This assumption is not restrictive, as better shown in example two, but, at first, it helps to keep the problem formulation simple and easy to manage. The priority number assigned to each load is always known and it changes on event occurrence. Finally, the supervisor task speed is at least ten times slower than the load's embedded control, if any.

MILP Problem Definition
Let us assume n as the number of loads under the control of the Supervisor, each one defined by the priority number ξ = (ξ 1 ξ 2 . . . ξ n ) ∈ N n positive integers, and P = (P 1 P 2 . . . P n ) ∈ R n + positive real numbers as the nominal power of each corresponding load. The Supervisor defines the load connections by avoiding the generator overload. In other terms, the Supervisor must guarantee that ∑ n i=1 P i ≤ P AV , with P AV the maximum available power assigned to the loads, calculated as (1), which is the difference between the nominal power (P nom ) and the power absorbed by loads (P u ) that are not under control.
The generator model can be expressed in simpler form by using the Thevenin's theorem [35]. At steady state, the model is further reduced to an ideal generator with in series a resistor. The steady-state generator output voltage is calculated as (2).
with being v g the open circuit voltage, R g the series resistance and i the generator current. The MIL-STD-704F defines the DC bus voltage limits [36]. By assuming as v min the minimum allowed voltage and i nom the maximum current corresponding to the minimum voltage, the generator nominal power for users can be defined by multiplying both members of (2) as: It must be noticed that the new nominal power definition comprises the maximum working power for the generator P nom < P gen−max . Being within the limits of (3) implies also the compliance of the minimum allowable voltage for the DC bus, thus giving to the MILP problem definition an homogeneous form, consequently overcoming its slack formulation.
The output is defined by the vector x = (x 1 x 2 . . . x n ) ∈ {0, 1} n corresponding to the status of the switches S = (S 1 S 2 . . . S n ), with the assumption of 0 for the open and 1 for the closed state.
The MILP problem is reported in terms of (4) with ξ and P defined above. Let us introduce a new priority vector ξ obtained by sorting the priority vector ξ in descending direction, with the assumption that a higher number corresponds to higher priority. The vector P is obtained from P with the condition of having the same index order of ξ. Thanks to the load properties and to the formulation of the problem as reported in (4), two considerations can be consequently stated: indicate the cumulative sum of the first P(k) elements, with k ∈ W as the integer positive number.

Corollary 1.
The cumulative sum of the vector power P increases monotonically.

Corollary 2.
The quadratic function (Σ P (k) − P AV ) 2 has only one minimum point k * .
The optimization problem (4) cannot yet be used to solve the above-defined problem illustrated in Figure 1, since some loads can be derated in power, while other ones allow only shedding. This issue can be solved with a different definition of the state variable x for loads allowing the power derating: in fact, by considering the load as a set of m small loads, each one having the same priority, the associated scalar state x i is substituted by the vector a vector having all elements equal to the scalar value associated to the element i. For most common problems, a division into m = 1024 set of power elements is acceptable for the application of our proposed methodology; in fact, by substituting the vectors into (4), we obtain a new set of variables: for example, if the third load is variable, the new vectors become: The power ratio for the load i is finally calculated as: Even if the problem dimension increases, however, the MILP output now directly gives the solution to the problem. As a drawback, the MILP solver has a critical issue due to the performance impact in terms of computational burden.
The optimization problem (4) can be formulated also as the knapsack problem and several algorithms are proposed in the literature to solve it [37]. In this paper, two different approaches are studied and compared to each other in terms of the steps needed to reach the optimal solution. It must be noticed that the number of steps is directly related to the performance in terms of computational time. The first considered method is the branchand-bound (B&B) one, commonly adopted to solve searching problems; it also allows the introduction of heuristics to speed up the convergence to the optimal solution. The second method uses the system knowledge to speed up the research of the sub-optimal solution and, finally, it refines the output by analyzing the slope error.

Branch-and-Bound Method
The branch-and-bound (B&B) algorithm is based on a tree search strategy and, in its original form, it divides the problem into a set of small sub-problems. This implies that the space solution is divided into smaller regions and the final solution is found recursively (i.e., branching); the rules established previously are used to prune off regions of the search space that are probably sub-optimal (i.e., bounding). For the problem (4), the binary branching formulation can be adopted.
The tree must be constructed respecting the priorities and the power level of each node, as shown in Figure 3a. In its more general formulation, the tree itself is constructed respecting the priorities of the nodes, and the search engine must calculate at runtime for each node whether the cumulative sum of the node power until the node under test satisfies (4). A small optimization is achieved by changing the top node, but the performance can be doubled at maximum. Therefore, a better approach must be adopted to reduce the computational time.
By assigning to each node the cumulative power Σ p (k), the tree can be constructed by considering the relationship between the cumulative powers. Then, (4) is implicitly defined with tree construction. Each node is associated with the subspace of the solutions; the left branch is related to the ones having higher cumulative power, the right branch to the smaller cumulative power. The search stops at a node having a greater value, respecting the condition expressed in (4). Despite the new solution, this still does not ensure that the tree is fully balanced, i.e., the left path has the same depth as the right path, and its structure is close to that of Figure 3b. The algorithm explores the tree until the cumulative power is greater than the available power P AV . As an example, we can refer to the to the algorithm of Figure 4a by assuming the tree nodes having as structure the left and right pointer, the cumulative sum value and the k index. The search time depends on the number of nodes that must be considered to find the optimal solution and finally depends on the tree balancing level.
As a case study, we consider the one reported in Figure 4b, having eight nodes and assuming P AV = 20 kW; starting from the top, the solution is searched in the right branch since the beginning nodes have cumulative power greater than P AV . Since the second node is smaller than the reference power, another test is made on its left node, but the test is not passed since it is associated with 25 kW as cumulative power. The optimal solution is given with only three steps, approximately equal to the depth of the position of the node d. Assuming that the tree is balanced, the number of steps required to reach the solution depends on the height of the tree or the depth of the k element: The worst case ε * is calculated as: with ceil being a function that rounds the number to the nearest major integer and log 2 x = log x/log 2.

Proposed Method-Linear Regression Approximation (LRA)
If all loads are equal, the Σ P (k) variable increases linearly and can be expressed with the form Σ P (k) = a · k. The minimum of the quadratic function: is k * and, in the case of all loads equal each other, is placed at the point k * 0 = P AV /a. The coefficient a can be approximated with a = Σ p (N)/N = P and it is equal to the average power of the loads; finally, the k value can be calculated as: If k * 0 = k * , the gradient method approach is adopted to find the solution. The algorithm searches the point that minimizes (8) and, thanks to the two corollaries previously reported, the local study of the function allows us to fully understand the side of the minimum with respect to the test point. To better understand this process, we can refer to the Figure 5, where the sign of the J variation is reported by changing the k index in both directions. By assuming as positive the direction associated with the k increasing and as negative the opposite one, starting from a generic point of the curve, the direction for the minimization of the function is given by the following relationship: The algorithm stops when J(k ± 1) > J(k). If the load's nominal power is relatively great compared to the available power, the solution is found in a few steps. For our study, then, we have to consider the worst case, i.e., several loads, with each one having small nominal power compared with P AV . In order to better understand this point, we first assume n = 80 small loads having power in the range [0 , 0.0467] P AV . Due to the unpredictable knowledge of the distribution of the power set, the Monte Carlo numerical approach is adopted to perform a statistical analysis. Figure 6 reports the stochastic test, repeated with 1000 different sets of vectors P, each one randomly generated. The J(k) function defined in (8) is shown in the Figure 6a subplot with the corresponding minimum point marked in red. The second subplot of the same figure reports the prediction error, calculated as: e k = k * − k * 0 ; the x-axis values are the corresponding error test numbers.
The mean error is found in one step and the worst case is in 9 steps. The worst case increases its steps when the loads have wide variation of nominal power; for this reason, a deeper study is needed to better understand this point.
Let us evaluate the prediction error of the random function e k (j) = k * (j) − k * 0 (j), with j being the power "granularity", i.e., the maximum power that each load can assume with respect to the available power, as defined in (11).
Each sample of Figure 7 is the statistical analysis of 1000 random numerical experiments, as conducted for Figure 6, by fixing the j value.
It can be observed that the mean error has a limited reliance on j and it is always less than 1.5. For j > 0.02, the variance and the maximum error also decrease, thus making the proposed method interesting as an alternative approach to the branch-and-bound one, or as a heuristic methodology for developing fast research algorithms.

Results
In this section, the branch-and-bound method is compared with the proposed linear regression one, in order to compare them in terms of time-computing performance, thus defining a new heuristic method to improve the Supervisor speed; the comparison results are reported in Table 2. Table 2. Methods comparison. d is the depth of the node, n is the number of elements, k 0 = 3.7, a = 0.19, j 0 = 0.017.

Method Mean Steps Maximum Steps
Branch-and-Bound d ceil(log 2 (n + 1)) As a conclusion from the comparison, we have that branch-and-bound method's step time is a function of the distance between the solution node and the top node. The proposed LRA method requires a mean step smaller than a constant, i.e., less than 1.5. By comparing both methods in the worst case, the branch-and-bound one needs an increase in the number of tests to reach the solution in a logarithmic relationship with the connected loads; in fact, tripling the steps, the number of loads must increase in the magnitude of one decade. Figure 8 reports in detail the comparison between the two methods. The graph is useful to select the best algorithm for the application by considering the worst case. Starting from the problem size n, the blue line gives the maximum number of steps needed to reach the problem solution. If we draw a horizontal line passing for this point, we obtain an intersection with the red line. The corresponding point gives the maximum allowable granularity j associated with the same maximum number of steps. If the problem under study has better granularity, the LRA method wins over the branch-and-bound method. As an example, if n = 10,000, the LRA method wins for j > 0.005. The LRA always wins in the case of j > 0.1 as its trend does not depend on the problem dimension, thus becoming simply a function of the granularity of the loads. The graphical method of Figure 8 can be also expressed in mathematical terms as follows: by comparing the maximum steps for the two methods, the condition making the LRA method perform better than the B&B one can be easily established, as defined in (12): For example, if j ≥ 0.02 p.u., the worst case is confined under 10 steps, as regards the branch-and-bound method, but with n * = 456, as shown following the gray dotted lines in Figure 8.
For those problems having a size of thousands of elements, the LRA method shows better performance with j > 0.02, since it monotonically decreases the number of steps by increasing j and it also shows the opposite convexity trend with respect to B&B.
The last aspect that we can consider for the comparison is the memory usage: the LRA algorithm needs a vector of n nodes, as the output is the vector index itself, while the B&B needs 4n nodes, i.e., three times greater than the LRA; for embedded applications, with a limited amount of memory, this is a critical aspect and it suggests that the LRS method should be preferred for all those applications having thousands of elements.
In summary, both the advantages and disadvantages of the two proposed methods are highlighted in Table 3. Table 3. Comparison between the two high-performance methods.

Method Strength Weakness
B&B computational time known a priori high memory usage time-consuming for lots of elements LRA low memory consumption variable execution time reduced computational time

Numerical Results
In this section, two numerical tests have been conducted, as better described in the following two subsections.
The generator feeds seven loads L1. . .L7, having the nominal power reported in Table 1, and external loads too, fully observable but not under the control of the supervisor. We refer to the one wire power flow circuit of Figure 9 and we assume P gen = P u + ∑ 7 i=1 P i for each timepoint. The power flow circuit is commonly used for studying large plants, also with three phase circuits; for these, each load is characterized by its active and reactive power [38]; for DC grids, instead, the loads are fully described by a scalar number representing their power.

L1 L2 L3 L4 L5 L6 L7
Other loads Generator P 1 P 2 P 3 P 4 P 5 P 6 P 7 P u P gen -ON/OFF loads -Programmable loads As for the test, we assume, as an example, that the loads L3, L4, L6 and L7 are variable with 128 steps. The remaining loads can be only switched OFF. With the proposed configuration, the problem dimension is n = (128 + 1) · 4 + 3 = 519. To reduce the postprocessing analysis, the load priorities increase from left to right, so ξ 1 > ξ 2 > . . . > ξ 7 . As better described in the second test, the assumption of knowing the load power, or of fixing it at the nominal value, is not restrictive, since the power variation can be considered outside the MILP formulation.

The Available Power Changes as Ramp
The first test has been conducted with the loads reported in Table 1, by assuming all of them at nominal power, as shown in Figure 9, in which the single-line diagram of the on-board plant is provided. Note that other loads (indicated with P u ) can be also connected to the plant, without any restrictions on the applicability of our proposed method.
The available power, defined as (1), changes as a ramp and it is shown by the blue line in Figure 10a. For the definition, the P AV is the maximum power that loads can absorb; in other words, the cumulative sum of load power, shown with the red line in Figure 10a, never can exceed the available power. Note that the power requested by the loads is always lower than the available power, both in ON/OFF and in the linear ones. The rules are always respected and the MILP always selects the configuration that maximizes the number of the supplied loads.  The ON/OFF loads are connected only in the case where the cumulative power is below the available power. The variable loads change their power according to a ramp trend. Figure 10b shows the load activation times and related p.u. power. By taking into account that the L1 load has the higher priority, while L7 has the lowest one, it can be observed that the activation sequence is always respected. It must be noticed that, for variable loads, the output of (5) has been considered instead of the binary load decomposition. As the test results conclude, the MILP formulation respects the available power limits, maximizes the absorbed power and respects the load priorities, thus feeding first higher-priority loads.

The Loads L1 and L2 Have Power Depending on the Time t
When one or more loads change their power, two approaches can be adopted. The simplest one is to update the cumulative power vector for each change, but, in this way, the solution is reached very slowly, thus reducing the system performance. A different approach, adopted in this subsection, generalizes both the load power and the available power.
Let us write the load power as: where P Lk is the constant term used for the power vector construction and p Lk (t) is the time-dependent term taking into account the power variation. By fixing the power vector P, the resulting available power can be calculated as reported in (14), i.e., as the difference between the available power (previously calculated using (1)) and the time-dependent term p LK . The new formulation is not limited to one load and it can be generalized for all loads only by changing the nominal power during the time of supervision. P AV (t) = P AV (t) − p Lk (t) (14) For the test, the loads L1 and L2 are assumed as time-dependent, as shown in Figure 11c. The blue line represents the instantaneous power of the loads, the activation time, the peak value of the requested power and, finally, the duration time depending on the loads; indeed, all these parameters are outside the control of the supervisor. For the MILP definition, the load nominal power has been considered, as shown with the red dotted line in the same figure; as a consequence, the p LK (t) can be calculated from (13) and it is shown as a yellow line in the figure. Finally, p LK (t) is used to define the available power P AV (t) in (14), as indicated by the blue line in Figure 11a.
With the new definition of the load power, the Supervisor allows the loads L1 and L2 to be connected, but this matter does not automatically imply that L1 or L2 are effectively absorbing power. For the algorithm test, it must be verified first that the absorbed power never exceeds the limits, which are now a function of the time-dependent load power. As can be observed, the absorbed power (red line) in Figure 11a is always less than or equal to the available power (blue line), so that the power limits are always respected. As regards the activation sequence, the available power is always greater than zero, thus allowing the first six loads to be always connected, as shown in Figure 11b. Only the low-priority load L7 has been used to avoid power overload by reducing its absorbed power.
(b) (c) (a) Figure 11. Test case for different available power and derating loads, with L1 and L2 changing with steps. (a) Available power compared with absorbed power, (b) p.u. setting for loads depending on properties: fixed loads and variable loads, (c) L1 and L2 power.

Discussion
This paper initially defines the problem of the power management for on-board electrical systems. In particular, our case study deals with a power management method installed aboard an aircraft. Our proposed method leads to the definition of the optimal configuration of the fed loads, so that both their priorities in supply needs (due to the flight necessities) and the generator power limit constraints are respected in any situation.
The proposed method's strength consists of describing the problem considering some loads as ON/OFF switches, i.e., they allow only shedding and some other variables.
Linear programming strictly defines a problem with integer value variables, but its adoption cannot be applied when they are both integer and continuous.
By dividing the large loads into smaller ones, each one described by only an integer variable (ON/OFF), it has been possible overcome to the limitation of power derating, but the growth of the problem dimensions has been consequently increased.
For all these reasons, two very high-performing methods have been compared, in order to evaluate the best one for the desired application (load feeding aboard an aircraft).
Many simulations have been carried out in order to validate the proposed method.

Conclusions
In this paper, the optimal power flow problem for a small priority-centered microgrid is proposed and solved by using the MILP theory.
Despite the fact that the MILP is promising for this class of problem, thanks to the reduced time needed to reach the optimal solution, its use for the DC power distribution systems is still a challenge due to the quadratic problem formulation when the grid is expressed in terms of the conductance matrix.
The power-balancing method for a steady-state grid is proposed, obtaining as a first benefit the linear formulation of the problem; moreover, the new model allows us to express the voltage limits in terms of power, giving a homogeneous formulation to the problem.
As an additional contribution, the MILP method has been extended also to variable loads; the proposed methodology allows us to manage both ON/OFF and variable loads using the same formulation; in fact, the variable loads are seen as tens of small ON/OFF loads and their contribution in terms of power is given by the number of small virtual loads connected to the bus over the size.
With this last assumption, the dimension problem consequently increases, so that the last contribution of the paper has necessarily to be the performance study of the two high-speed algorithms here presented. Finally, the paper considers loads having an instantaneous power function of the time; it has been reported an additional method that does not need any change in the cumulative power vector, thus keeping its performance still satisfactory in terms of computational burden.
Future works will be related to the study and the implementation of real-time performance tests using different-order smart grids, taking into account also the transient phenomena.
Author Contributions: L.R. developed the theoretical formalism and performed the analytic calculations. G.R. and P.C. performed the numerical simulations. All authors equally contributed to the final version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: