MPC Framework for the Energy Management of Hybrid Ships with an Energy Storage System

: This paper proposes an advanced shipboard energy management strategy (EMS) based on model predictive control (MPC). This EMS aims to reduce mission-scale fuel consumption of ship hybrid power plants, taking into account constraints introduced by the shipboard battery system. Such constraints are present due to the boundaries on the battery capacity and state of charge ( SoC ) values, aiming to ensure safe seagoing operation and long-lasting battery life. The proposed EMS can be used earlier in the propulsion design process and requires no tuning of parameters for a speciﬁc operating proﬁle. The novelties of the study reside in (i) studying the impact of mission-scale effects and integral constraints on optimal fuel consumption and controller robustness, (ii) benchmarking the performance of the proposed MPC framework. A case study carried out on a naval vessel demonstrates near-optimal and robust behaviour of the controller for several loading sequences. The application of the proposed MPC framework can lead to up to 3.5% consumption reduction due to utilisation of long term information, considering speciﬁc loading sequences and charge depleting (CD) battery operation.


Regulatory and Technical Context for Battery Installations
Reduction of the environmental impact of ships has been a subject of concern and has been addressed by regional and global regulators. According to the fourth International Maritime Organisation (IMO) greenhouse gas (GHG) study in 2020 [1], shipping contributes an average of 2% of total anthropogenic CO 2 emissions, while with business as usual (BAU) projections, the CO 2 emissions are expected to increase by 0-50% compared to the 2018 figure. To that end, IMO has been actively involved in, amongst other activities, enhancing the energy efficiency of ships, within the MARPOL Annex VI regulatory framework. The IMO's efforts aim to reduce GHG emissions to 50% of the 2008 baseline value by 2050. According to all future projections in the study, it is difficult to achieve the CO 2 reduction goal solely by means of more efficient conventional propulsion and lower sailing speeds. Towards that end, several alternatives to conventional shipboard energy generation have been proposed in the literature. Listing only some of the proposed alternatives, reduction of emissions compared to conventional power plants has been reported by using liquefied natural gas (LNG) fuel [2,3], ammonia fuel and ammonia injection [4,5], biofuels [6,7], and energy storage systems (ESSs-list of abbreviations given in Table A2).
Although ultracapacitors are utilised when surges of power are needed by electrical consumers on-board (e.g., weapon systems on naval vessels) [8], the scope of this study is limited to batteries, as the most researched and most promising ESS in the maritime sector [9,10]. The price of lithium-ion batteries has plummeted over the past 10 years, by 1.
The shape of the efficiency curve of the prime mover is the dominant decisive factor for the applicability of an ESS installation 2.
When an ESS replaces engine power instead of adding power to an existing engine rating, the efficiency gains are smaller especially at nominal loads 3.
In electric propulsion, most conversion losses introduced by an ESS are already in place 4.
A DC grid allows for use of non-fixed-speed generator sets and consequently affects the powertrain efficiency in favour of using generator sets and not ESSs 5.
Efficiency gains in electric propulsion due to ESSs are found to be 4.7-7.8% at 50% load and 18-30% at 15% load.
Hybrid propulsion is also found to yield efficiency benefits when combined with a battery unit in charge depleting (CD) mode (shore charging) [16]. It has had commercial success in tugboats by Damen and yachts by Feadship [9]. In [17], four different propulsion configurations in megayachts have been compared and significant advantages in hybrid propulsion have been identified when combined with fast shore charging in marinas with local smart grids. The feasibility of zero-emission high-speed catamaran vessels has been one of the subjects for [18,19], where the importance of charging stations has also been identified as crucial. Finally, battery installations can be combined with fuel cell systems on board ships in order to improve efficiency and dynamic behavior [20,21].
The environmental analysis of battery installations should not be limited to efficiency benefits. An approach taking into account different impact categories, accurate measurements and the life-cycle (life cycle assessment) of the vessel is capable of providing full ESS installation [10] it is potentially beneficial for an EMS to be deployed in low fidelity systems during early design.
In model predictive control (MPC) the future disturbance estimation over a set prediction horizon span is made online, then a quadratic programming (QP),dynamic programming (DP), or other optimisation solution is used to find the optimal control sequence for a control horizon span, which may not necessarily be of the same length [23,29,30]. MPC is computationally expensive and requires an accurate model description of the system to yield satisfactory results, while, in order to be effective, prediction techniques must be used [31]. In [32], an MPC EMS is proposed, utilising intelligent transportation systems (ITS) for driving cycle prediction. A combination of route prediction from the telematics system together with MPC for powertrain control has also been implemented in [33]. Online implementations of predictive control with information supplied by the GPS and the traffic-flow information systems yield near-optimal results in [34].
Shipboard MPC implementations may function as either power or energy management strategies depending on their timescale and objective function. If the timescale of the control horizon is in the order of several seconds and total fuel consumed is not in the objective function, the strategies can be classified as power management strategies. Such controllers take into account the transient behaviour of the powertrain and ship system. For instance, in [35], the power demand fluctuations (in the order of few seconds) in various sea states due to the propeller load are mitigated, and in [36], the nonlinear ship response of an offshore support vessel due to environmental disturbances is taken into account to control ship speed in the surge direction. Another power management implementation is found in [37].
In contrast to power management strategies, there are only a few marine studies on shipboard energy management for entire missions. In [26], a predictive controller that uses a short-term prediction scheme has been implemented in an electric tugboat, yielding 9% fuel consumption reduction when compared to a rule-based EMS. The study also incorporates the SoC constraints in the optimisation horizon solver. However, mission-scale predictions surpassing the prediction horizon and their effect on optimal fuel consumption were not examined. Furthermore, a comparison with optimal control policies was not present. Huotari et al. [38] takes into account mission-scale predictions and benchmarks the proposed MPC-based controller against an exhaustive solution for a diesel-electric cruise ship power plant case study. This study uses mixed integer linear programming (MILP) as the MPC solver while it also incorporates a 2-stage predictive scheme to provide a mission-scale prediction. Despite the effectiveness of the method in achieving close to optimal results, it is noted that the controller does not incorporate updates on the future disturbances while sailing and the prediction is provided to the solver prior to the mission. A rule-based controller is utilised instead in order to compensate for deviations between the actual disturbances and the initial estimate. In addition, a constraint on the maximum and minimum (dis)charging power of the battery is incorporated in the optimisation problem. This should be effective in protecting the battery from exceeding the minimum SoC value. However, a constraint on the SoC value at the end of the mission is not present, which is important for optimal performance in charge-sustaining mode.
In an automotive study [34], an implementable MPC framework is developed with near-optimal performance, enabling online updates and utilising information from the ITS for load prediction. This controller is adapted in the present work to design a shipboard MPC framework and address the performance issues of advanced EMS in ships.

Aim, Contributions and Assumptions
The first goal of this research work is to examine how the loading sequence of the whole mission affects the fuel savings of an exhaustive solution found with dynamic programming. This is attained by using partial and full loading sequence information, corresponding to short-term and mission-scale predictions made for the EMS.
Secondly, this study develops an energy management framework based on MPC that addresses multiple objectives which have not simultaneously been fulfilled by EMS alternatives proposed in the literature. The MPC framework:

1.
Deals with battery-induced integral state constraints robustly.

2.
Realises potential fuel savings due to mission-scale information.

3.
Performs close to the exhaustive DP solution results.

4.
Can be readily set to function in charge-sustaining (CS) or charge-depleting (CD) battery mode.

5.
Enables real-time adaptation to shipboard updates on the mission-scale disturbance estimation while sailing, rather than using offline tuning of parameters such as the equivalency factor in ECMS, predefined power demand sequences [38] or rule-based tuning. 6.
Requires minimal parameter tuning on specific operating profiles enabling deployment on early design stages, with the prerequisite that a model description for the powertrain is available.
This MPC framework has been validated in both CS and CD mode by incorporating the controller in a dynamic Simulink ® model of a naval hybrid propulsion powertrain with hybrid power supply, validated in [39].
It is noted that, to factor out the performance of individual long-term and short-term prediction solutions, the study assumes perfect predictors with a fabricated information barrier scheme for the experiments, which will be discussed in more detail in the results segment. For future reference, several alternatives for disturbance prediction found in recent studies are proposed. Long-term predictions can be made by data-driven prediction approaches, such as the one found in [40], where a physics-based machine learning (PBML) model is proposed. The utilization of advanced routing models, such as the one found in [41], can be combined with weather data, navigation aids, the governor's human input, and sensor information to provide highly accurate predictions. For short-term disturbance estimation, several prediction schemes may be used, such as the prediction scheme proposed in [26], neural network approaches, such as the RBF-NN implementation in [34] and autoregressive moving-average (ARMA) models [42].

Model Description
The propulsive plant is described by a set of differential and algebraic equations (DAE) with inequality constraints for the operating limits of the various components. The DAE system is defined by the propulsive layout and by means of fitted equations whose parameters are tuned using experiments on a component level as in [16], either directly from physical measurements or from detailed dynamic simulation models of components. For fuel consumption minimisation, quasi-static models are typically modelled by fuel consumption maps, or fitted functions, as dynamic behaviour is of secondary importance [30]. Most of the formulations for the components have been verified and validated in [16] and the derivations will not be repeated here.
The quasi-static description for the components of a shipboard hybrid power plant can in turn be reused for different hybrid layouts. Following the methodology for formulating a non-linear model description for MPC purposes as in [34], the model description for the hybrid propulsion hybrid power supply configuration in Figure 1 is presented here: 1.
The electrical and mechanical nodes of the plant correspond to energy balance equations.

2.
For a given disturbance vector, which is typically comprised by the propulsive load, hotel load and shaft speed, the operating setpoints of the components are not uniquely defined due to the degrees of freedom (DoFs) in the system. For example, the propulsive demand can be met by infinite combinations of main engine and induction motor power output levels. The number of indefinite variables in the model description corresponds to the DoF in the system and each DoF correspond to an additional system control input.

3.
Each component introduces operating constraints which should be explicitly expressed.

4.
Finally, since the SoC of the battery must be within its operating constraints, it should be introduced to the system as a state variable.
1. The electrical and mechanical nodes of the plant correspond to energy balance equations. 2. For a given disturbance vector, which is typically comprised by the propulsive load, hotel load and shaft speed, the operating setpoints of the components are not uniquely defined due to the degrees of freedom (DoFs) in the system. For example, the propulsive demand can be met by infinite combinations of main engine and induction motor power output levels. The number of indefinite variables in the model description corresponds to the DoF in the system and each DoF correspond to an additional system control input. 3. Each component introduces operating constraints which should be explicitly expressed. 4. Finally, since the SoC of the battery must be within its operating constraints, it should be introduced to the system as a state variable.
and by discretising with a timestep of Δt: where Icell,bat the current in the battery pack (positive for discharging). The SoC of a battery cell is defined as the charge in the battery, divided by the cell's maximum charge,Q(t)/Q cell,0 . The derivative yields and by discretising with a timestep of ∆t: where I cell,bat the current in the battery pack (positive for discharging).
There are several ways to define battery efficiency, categorised in [43] into global and local efficiencies, depending on whether the efficiency is averaged over the entire (dis)charging cycle (energy ratio) or is calculated as a power ratio. Local efficiency varies by both (dis)charging battery power and SoC. Local efficiency is strongly dependent on the (dis)charging battery power, while dependency on the SoC is much weaker [43]. This is the reason why global efficiency is sufficient for CS operation. However, in CD battery mode, the SoC value is depleted at 20% [16], and in such low values contribution to the battery efficiency is no longer negligible. The battery efficiency is used in EMSs for the estimation of the SoC and small errors are accumulated over time.
This study uses the battery equivalent circuit found in [16]. However, instead of the global Ragone efficiency, for the reasons above, the electrochemical local efficiency of the battery η ech has been used instead. The electrochemical power of the battery cell is defined as the power of the battery without taking into account the internal losses [43]: where V cell,open the open circuit voltage. The electrochemical power P cell,ech is connected to the battery cell power P cell,bat with the electrochemical local efficiency: The electrochemical efficiency can be expressed as a polynomial function of the power and SoC: where c ij and d ij are fitting parameters for charging and discharging to be determined experimentally, P n = P cell /P cell,max . A complex expression for the open circuit voltage is used as in [16]: where α i model parameters to be determined experimentally. Combining Equations (3) and (4) while also converting from the pack to the cell: where n par is the number of cells connected in parallel and n ser in series. The SoC for every the next time-step now is: Substituting the expression for the I cell,bat into Equation (8) provides a way to calculate the SoC for the next time step.
For non-DC architectures, the battery pack is connected via a static converter (inverter/rectifier) to an AC grid. The battery module current I bat,AC , for the case of a3-phase AC grid, is connected to the battery pack power P bat,AC with the expression: where V line is the line voltage and pf is the power factor of the AC bus. The associated current losses at the converter, I bat,cnv,los , are given by the following quadratic fitted relationship, assumed to stay the same for both converting modes: I bat,cnv,los = I nom,los · e 1 I 2 n + e 2 |I n | + e 3 , where I nom,los the nominal current losses and e i fitting parameters. Moreover, I n is the normalised current I n = I bat,AC /I bat,AC,nom : The I pack,bat can then be found as: I pack,bat = I bat,AC + I bat,cnv,los .
Note that I bat,cnv,los > 0 so that during charging (I < 0) the rectifier losses are subtracted from the current input towards the battery, while during discharging (I > 0) the inverter losses are added.
The fuel consumption of the diesel generator set can be measured by varying the torque at a constant rotational speed. A cubic approximation is provided using manufacturer's data (as in [16]) as a function of the normalised power output of the generator set, P genset,n : .
where g i positive fitting parameters and k genset the number of active generator sets at a time. The approximation is a convex function because the second derivative of the fuel consumption is 2g 2 > 0, strictly positive. This allows for choosing the optimal number of active generator sets using only the required power output (See Appendix A). The main diesel engine fuel consumption can be approximated with a polynomial function of the normalised power output P DE,n and engine shaft speed n DE,n : .
where b i,j fitting parameters, k p the number of propulsion trains (assuming one engine per propeller). The electric machine has been modelled based on a piecewise polynomial fit, where the normalised power losses, P PTO/I,loss,n = P PTO/I,loss /P PTI,loss,nom , of the variable speed drive can be expressed as a symmetric function of the normalised motor speed, n n , and torque, M n : where v i parameters. The power of the electric machine is P PTO/I = (T nom N PTO/I,nom )2πM n n PTO/I,n , with T nom and N nom as the nominal torque and shaft speed, correspondingly, and n PTO/I,n the normalised PTO/I shaft speed. If for power-take-off, P PTO/I < 0, and for power-take-in, P PTO/I ≥ 0, the power before the PTO/I gearbox is: while the power at the AC grid is: P PTO/I,AC = −P PTO/I + P PTO/I,loss .
The gearbox ratio between the propeller shaft and the PTO/I, and the main engine and the propeller shaft is: The mechanical power balance can be expressed as where P prop is the total propulsive demand and η GB is the gearbox efficiencies for the PTO/I and the main engine shaft. The power quantities have been multiplied with their nominal values. The power balance at the AC grid is k genset P genset + k p P PTO/I,AC + P bat,AC = P hotel , where P hotel is the hotel power demand on the AC grid and the k genset the active generator sets at each instant.
Since the model is quasi-static, the dynamic behaviour of the propulsion system due to waves, propeller shaft inertiae, and propeller pitch are neglected. Thus, the propulsive power demand P prop and the propeller shaft speed n p can be separate disturbances to the system. A separate module that generates these disturbances from e.g., the sea state and the ship speed can be used.
The operating envelope of the main engine is divided into three different regions, the power limit, n 2 ≤ n DE ≤ n max , the torque limit, n 1 ≤ n DE ≤ n 2, and the turbocharger limit, n min ≤ n DE ≤ n 1 . The last limit is expressed in terms of maximum power as a function of the shaft speed: where β 1 , β 2 parameters. A safe factor can be used in the controller to ensure that the operating limits are not violated during operating point transitions, and the engine is stable.
The per cell maximum continuous current of the battery is I max , while the cell capacity is Q 0,cell . This gives a maximum C-rate of C max = I max /Q 0,cell per cell. For this C-rate, the power output P cell,max can be measured, resulting in a total power output of P bat,max = n ser n par P bat,max . Furthermore, there is a limit on the battery state of charge, SoC min ≤ SoC ≤ 1. It is important to note that the latter is an integral constraint on the state.
A minimum value for the torque M n,min of the electric machine is active, or zero when the electric machine is disengaged.
Finally, the power from each generator set is P genset,min ≤ P genset ≤ P genset,max or zero when the generator set is off.

Energy Management Strategy Framework
The ship dynamic behaviour can easily be written with the system description of the following form: S where SoC is the state variable, y the system output vector, u the control input vector and r the disturbance vector, which are: The algebraic power balance in Equations (21) and (22) allow for the system output to be calculated. Note that the physical description contains non-linear expressions, e.g., due to the (dis) charging of the batteries. The operating constraints presented previously are part of the system description.
A non-linear constrained optimal control problem (OCP) is introduced where the objective is to find the control input at every time instant which minimises the cost function J: where the fuel consumption terms and the G term are functions of the control input vector, the state, and the time. Term ϕ penalises the end state SoC and term G penalises the state variable. t f and t 0 the end and start instants. Discretising with respect to time, the OCP then requires a control input policy π* = {u 1, u 2, . . . , u N },u k ∈ U so that: where the constraints on SoC are enforced with term γ k . Time steps t = t 0 + 0∆t, t 0 + 1∆t, . . . , t 0 + N∆t are denoted as 1,2, . . . N. Every control action results in the state of the next time step. Further discretising the state and control input variables into N SoC and N u values, correspondingly, allows for the OCP to be solved with dynamic programming DP. Using the MATLAB ® generic DP solver in [44], the method has been followed in a number of studies for automotive applications to minimise J [30,34,45]. The momentary fuel consumption for the generator set and the main engine is given in Equations (13) and (14), correspondingly, and is present in the objective function of the OCP in Equation (31). The genset fuel consumption is only dependent on the demanded power from the generator P genset (see Appendix A), which is in the algebraic Equation (22). The main engine fuel consumption depends on the demanded main engine power, P DE, and the shaft speed n s , which is one of the three disturbances to the system. For a valid solution to the model description with given disturbances, the value for the PTO/I power is found from Equation (16) for torque control input, M n , and the battery power is found from Equation (9) for battery current control input, I n, are found. Together with the disturbances for the hotel, P hotel , and propulsive power demand, P prop , it yields two Equations (21) and (22) with two unknowns: the power of the main engine and the power of the generator set. Then from Equations (13) and (14) the solution can be evaluated with respect to the objective function in Equation (31). The solver finds the values for the control inputs which minimise the momentary fuel consumptions, summed up for all time steps.
Operating limits on the components of the ship's hybrid powerplant and the integral constraints on the SoC need to be applied effectively. The DP solver is using only the defined set of the discretised control input vector U k . It also includes the penalising terms ϕ and γ to apply constraints on the state. This allows for robust behaviour, ensuring all constraints are enforced.
In -MPC-, a series of OCPs are solved for subsequent time intervals or control horizons T c = N c ∆t. At each time step k,where t k = k∆t the state variable x k is sampled from the physical system. The disturbance vector r k,k + Nc for the horizon needs to be known or estimated in advance. Each OCP solution is found by feeding the above to the system description, yielding the optimal control input policy for the horizon π * k,k + Nc = {u k , u k+1 , . . . , u k + Nc }. However, only the first control input vector u k is given to the physical system, while the rest are discarded and the procedure is repeated for the next time step k + 1, sampling the next state variable x k+1 from the system. The advantages of MPC are identified in [29]. MPC is capable of using multiple variables for the state and control inputs; it is also suitable for handling binary variables and non-linear systems. For non-linear systems, the term non-linear MPC is used in the literature [46]. Applicability is promising for shipboard hybrid systems, where more components are present compared to automotive hybrid powertrains, and (dis)connecting generator sets or (dis)charging the battery module introduces non-linear equations. Moreover, MPC is inherently able to use new predictions for the disturbance vector r whenever available, because it fully recalculates the OCP at every time step.
Conventional MPC originates from the process industry where it is important to use a setpoint tracking term to penalise deviations from the reference system output or state at each stage of the OCP [29]. However, the objective of an EMS fits the description of economic or performance optimisation and can therefore be addressed by formulations described in the literature as economic MPC [46,47]. These formulations drop the conventional setpoint tracking term and instead optimise for an economic objective by use of a performance term (first term of Equation (29)) [48,49]. The setpoint tracking is then done with a terminal constraint or penalty (second term of Equation (29)). The reader is referred to [46] for an overview on the stability and optimality of economic non-linear MPC.
Sun et al. [34] in a study for automotive vehicles, use a non-linear MPC approach as a compromise between DP and ECMS for energy management. On one hand, DP requires the load profile and duration to be known in advance, thus, it is predominantly used for benchmarking or tuning other EMS. On the other hand, ECMS performs well when information is given only on the present time instant, however, tuning of the equivalence factor, which depends on the loading profile, is not a trivial tasks and needs to be done in advance. This introduces a loading-sequence-specific control design overhead. Adaptive implementations (A-ECMS) use only a set of equivalence factors and the implementation in [16] failed to perform well for loads that differ for the profiles used for tuning.
Following the framework design found in [34], the proposed controller MPC framework is comprised by two modules, namely the reference trajectory generation (RTG) module and the non-linear MPC module (top two blue blocks in Figure 2). Both modules solve a discrete-time finite horizon non-linear OCP with DP in order to minimise the fuel consumed in the form of an integral objective function, given by Equations (30) and (31). In predefined time intervals, the RTG module is initiated and fed with the available remaining mission load estimation as the disturbance vector input r mission . This future update is shown in Figure 2 with the arrow from the mission-scale disturbance estimation at the top left to the RTG module at the top right. Every time a new future update on the remaining mission disturbance vector is fed to the framework, this also updates the RTG solution. From the RTG solution, a mission reference state trajectory for the battery SoC ref (t) is computed. The MPC module is instead using a short term prediction r horizon and is solving a new OPC for a defined control horizon at every control time-step t k . Using an end-state constraint (ESC) for the SoC at the end of each horizon, the MPC solution is guided towards the solution generated by the RTG module: where E is the parameter in the ESC tuned to allow for feasible solutions in case of disagreement between the short-term MPC solution and the mission-scale solution by the RTG module. The end-state constraint is enforced with the ϕ term of Equation (31).
(e.g., captain's input, weather and routing data), can be combined with a short-term disturbance estimation scheme, examples of which can be found in [26,34,42,50,51]. The use of the quasi-static model in the MPC module provides past values of the state and the system outputs that may enhance the performance of the estimation scheme. The present study uses a perfect prediction scheme for both the mission-scale and the short-term, so that the performance of the scheme does not interfere with benchmarking the performance of the proposed MPC framework.  In both the RTG and MPC modules the behaviour of the propulsive plant is described by the quasi-static model description given at the beginning of this chapter.
By using two different time-scales for the estimations, mission-scale information (e.g., captain's input, weather and routing data), can be combined with a short-term disturbance estimation scheme, examples of which can be found in [26,34,42,50,51]. The use of the quasi-static model in the MPC module provides past values of the state and the system outputs that may enhance the performance of the estimation scheme. The present study uses a perfect prediction scheme for both the mission-scale and the short-term, so that the performance of the scheme does not interfere with benchmarking the performance of the proposed MPC framework.

Model Description Tuning and Verification
The MPC framework was validated by incorporating the controller in a dynamic Simulink ® model of a hybrid power plant of a naval vessel provided by Damen Naval BV ( Figure 1). For that reason, the parameters of the system description have been formulated using measurements from the model, whose details and validation are mostly covered [43], except from the electric machine, which has been simplified here to the same polynomial in Equation (15).
The telegraph position, Tg, connects to the virtual shaft speed by the equation N virt = (N virt,max − N virt,min )Tg + N virt,min , where N virt,min and N virt,max are the min and max virtual shaft speed limits. The virtual speed is connected to the propulsive load and actual shaft speed for a fixed combinator curve. An interpolated fit has been chosen, shown in Figure 3. It has been used to generate the disturbance for the shaft speed and the propulsive load from the telegraph position. same polynomial in Equation (15).
The telegraph position, Tg, connects to the virtual shaft speed by the equation Nvirt = (Nvirt,max − Nvirt,min)Tg + Nvirt,min, where Nvirt,min and Nvirt,max are the min and max virtual shaft speed limits. The virtual speed is connected to the propulsive load and actual shaft speed for a fixed combinator curve. An interpolated fit has been chosen, shown in Figure 3. It has been used to generate the disturbance for the shaft speed and the propulsive load from the telegraph position. In Table A1the parameters of the vessel's power plant are given. On a component level, Figures A1-A3 present the fitted results for the electrochemical efficiency of Equation (5) and the main engine consumption of Equation (14). The Simulink ® model uses a simple polynomial expression for fuel consumption of the generator set as in Equation (13) so the reduction is trivial. The layout of the hybrid propulsion hybrid power supply power plant is given in Figure 2.
For the electrochemical efficiency of the battery, the effect of the state of charge is observed in Figures A1 and A2 but is more pronounced in the case of cell discharging. The biggest fitting error is observed in Figure A2 with a maximum value of 0.7%, at the region where the state of charge is above 0.95 and for a discharging power above 40%. The residuals for charging are shown in Figure A1, with the biggest fitting error of 0.7% occurring at the zone where the charging power is above 90%.
The bottom plot of Figure A3 includes the envelope boundaries and shows that the maximum residuals are observed alongside the operating limits of the engine, but do not exceed the 0.008 mark.
The fitted equations for the main engine (Equation (14)) and the battery cells (Equation (5)) are normalised in order to ensure that the fit is as accurate as possible. The corresponding nominal values for the normalised quantities are found in Table A1.
The model description has been verified offline by plotting the results of the same loading profile of total duration of 3000 s against their Simulink ® model counterparts, as shown in Figure 4. For positive values of the battery power in Figure 4b the charge of the battery is depleted, shown in Figure 4f, as described in Equations (6)- (8). For this inte- In Table A1 the parameters of the vessel's power plant are given. On a component level, Figures A1-A3 present the fitted results for the electrochemical efficiency of Equation (5) and the main engine consumption of Equation (14). The Simulink ® model uses a simple polynomial expression for fuel consumption of the generator set as in Equation (13) so the reduction is trivial. The layout of the hybrid propulsion hybrid power supply power plant is given in Figure 2.
For the electrochemical efficiency of the battery, the effect of the state of charge is observed in Figures A1 and A2 but is more pronounced in the case of cell discharging. The biggest fitting error is observed in Figure A2 with a maximum value of 0.7%, at the region where the state of charge is above 0.95 and for a discharging power above 40%. The residuals for charging are shown in Figure A1, with the biggest fitting error of 0.7% occurring at the zone where the charging power is above 90%.
The bottom plot of Figure A3 includes the envelope boundaries and shows that the maximum residuals are observed alongside the operating limits of the engine, but do not exceed the 0.008 mark.
The fitted equations for the main engine (Equation (14)) and the battery cells (Equation (5)) are normalised in order to ensure that the fit is as accurate as possible. The corresponding nominal values for the normalised quantities are found in Table A1.
The model description has been verified offline by plotting the results of the same loading profile of total duration of 3000 s against their Simulink ® model counterparts, as shown in Figure 4. For positive values of the battery power in Figure 4b the charge of the battery is depleted, shown in Figure 4f, as described in Equations (6)- (8). For this integrated verification case the setpoints do not yield an extensive charging period for the battery. However, the battery charging mode is verified also on a component basis, an important step for accurate SoC tracking, as it is later shown in Section 3.4 that the proposed framework yields extensive charging periods (see Figure 8d). The offline and online behaviour have their biggest difference of −4.1% for the diesel generator power around the 800 s mark and for a relatively short duration of 200 s. Apart from this region, the difference is well below 1% and it is concluded that the model description is accurate enough to replicate the Simulink ® model. It is clarified here that in this section only a model verification is presented and no controller has been implemented in the Simulink ® model. That is, the same control input and disturbance vector sequences were fed to both the model and its reduced quasi-static model description, comprised of all components described in the beginning of this section. around the 800 s mark and for a relatively short duration of 200 s. Apart from this region, the difference is well below 1% and it is concluded that the model description is accurate enough to replicate the Simulink ® model. It is clarified here that in this section only a model verification is presented and no controller has been implemented in the Simulink ® model. That is, the same control input and disturbance vector sequences were fed to both the model and its reduced quasi-static model description, comprised of all components described in the beginning of this section.

Verification of the Reference Trajectory Generator Module
The MPC framework was first implemented offline in order to select suitable values for the solver parameters and the MPC. Both the RTG and MPC modules of the framework use the generic MATLAB ® DP solver [45] to solve the OCPs.
Regarding the numerical details, the discretisation of the state and control input variables was found to be adequate at N x = N u = 41, above which not much reduction in fuel consumption was observed, while the computational load increased considerably. The number needs to be odd so that the mid setpoint of each variable can be zero, while the variables are normalised. The solver's boundary line method was used, which provides better accuracy for the discretisation number. The non-fixed grid option was selected and the time step of the solver was set to 20 s.
The RTG module underwent several qualitative verification tests: A base case of the telegraph position, an increased position by 10%, and a decreased position by 10% were compared to check whether the fuel consumption changed consistently. The results are shown in Table 1. The telegraph position is connected linearly to the propeller speed (until 0.9n p,max ) and the propeller speed is connected almost linearly to the propulsive demand (until 0.68P prop,max ) due to the controllable pitch propeller (CPP), as shown in Figure 3. The specific fuel oil consumption is fairly uniform (less than 7% difference between minimum and maximum value in Figure 5a). Thus, an increase in fuel consumption is expected with 10% increase in the telegraph position and vice versa, since that corresponds to more than a 20% increase/decrease for the power demand.

Verification of the Reference Trajectory Generator Module
The MPC framework was first implemented offline in order to select suitable values for the solver parameters and the MPC. Both the RTG and MPC modules of the framework use the generic MATLAB ® DP solver [45] to solve the OCPs.
Regarding the numerical details, the discretisation of the state and control input variables was found to be adequate at Nx = Nu = 41, above which not much reduction in fuel consumption was observed, while the computational load increased considerably. The number needs to be odd so that the mid setpoint of each variable can be zero, while the variables are normalised. The solver's boundary line method was used, which provides better accuracy for the discretisation number. The non-fixed grid option was selected and the time step of the solver was set to 20 s.
The RTG module underwent several qualitative verification tests: A base case of the telegraph position, an increased position by 10%, and a decreased position by 10% were compared to check whether the fuel consumption changed consistently. The results are shown in Table 1. The telegraph position is connected linearly to the propeller speed (until 0.9np,max) and the propeller speed is connected almost linearly to the propulsive demand (until 0.68Pprop,max) due to the controllable pitch propeller (CPP), as shown in Figure 3. The specific fuel oil consumption is fairly uniform (less than 7% difference between minimum and maximum value in Figure 5a). Thus, an increase in fuel consumption is expected with 10% increase in the telegraph position and vice versa, since that corresponds to more than a 20% increase/decrease for the power demand.   The second test introduces narrower operating boundaries for the components of the power plant. The minimum engine power has been moved from 45% of the nominal value to 50%, resulting in a consistent increase in fuel consumption. Likewise, an increase in the minimum torque for the PTO/I also increases fuel consumption. Figure 5 includes a series of plots that show how the various quantities of the components evolve for the increased minimum engine power value. It is noted that not all of the battery charge is used (as opposed to the base case) because the minimum main engine power plus the minimum PTI power is higher than the propulsive demand, meaning that only in peak propulsive power demand regions the PTI can be activated to further deplete the battery.

Verification of the Model Predictive Control Module
At the end of Section 2.2 an overview of MPC is provided. The implementation of the MPC module for this case study is solving an OCP for subsequent control horizons using a DP solver with an ESC at the end of each horizon. In more detail: 1.
The state is the battery SoC k and it is sampled from the Simulink ® model of the propulsive plant at each time step t k and used as the initial SoC for the OCP.

2.
The disturbances to the system are the shaft speed, propulsive power demand, and hotel power demand. The predicted values for each disturbance and for the scope of each horizon are provided to the OCP at each time step t k .

3.
The SoC MPC at the end of the horizon t k + N c ∆t is constrained by the SoC ref of the global solution at the same time-step. This is the ESC of Equation (32).

4.
The DP solver is utilised to provide the solution to the OCP above for the control horizon. The control inputs to the system are the normalized battery current I n and electric machine torque M n . The resulting control input vector for the next time step u k + 1 is fed to the system and the control input vectors for the rest of the time-steps, u k + 2 . . . Nc , are discarded.

5.
At the next time-step the MPC module samples the next SoC and iterates the procedure.
Before implementing the MPC module in the Simulink ® model of the propulsive power plant, the parameters in the MPC code are tuned offline so that the generated control inputs agree with the control inputs of the global solution, which is the DP solution for the full mission, provided online by the RTG module. In this offline verification, the SoC reads the initial state from the results of the global solution at each time-step. The goal is to ensure that the OCP solution with an ESC at the end of the control horizon yields the same control inputs with the global solution.
The ESC parameter for the MPC module was found offline by decreasing the value until the load profile provokes an infeasibility error. For E = 0.01 the framework is stable with 41 grid elements (N u = N x = 41), while E = 0.07 can be used for 61 grid elements, and E = 0.05 for 81. An analysis of the computational cost of the DP solver used by the MPC framework in this study can be found in [44]. The MPC module control horizon was found offline by starting at the short horizon of only two time steps (40 s) and increasing until the control inputs from the MPC module's OCP converge to those from the RTG module. It was found that a horizon of 600 s yields optimal response with no instabilities for both step-up and step-down load sequences ( Figure 6). The prediction horizon has been chosen to be the same as the control horizon.
for both step-up and step-down load sequences ( Figure 6). The prediction horizon has been chosen to be the same as the control horizon.

Validation of the Model Predictive Control Framework
The validation of the MPC framework is done by incorporating the controller in the Simulink ® model, using lvl-2 custom S-functions. A series of realistic telegraph trajectories, generated by a probability density function from historical data, was provided by Damen Naval BV for the vessel. Using the combinator curve presented in Section 3.1, these telegraph trajectories correspond to the disturbances for the propulsive load and the shaft speed. The hotel load power demand has been constant at 1000 kW.
The above loading sequences were spliced together to construct two different scenarios, referred to as Profiles 1-2. The profiles are used in an information barrier scheme which has been devised for the validation of the MPC framework: Two online optimal DP controllers are being fed the same disturbances with the MPC framework controller. At a certain time-step, called the information barrier point (IBP), the telegraph trajectory changes from the initial prediction to the actual trajectory (solid and dashed lines, respectively, shown for profiles 1 and 2 in Figure 7):

•
The first online DP controller has perfect information about the disturbances, including the change at the IBP. This means that it is the optimal controller and calculates the control inputs for the telegraph trajectory spliced at the barrier point.

•
The second online DP controller solves for the initial prediction until the IBP, where it recalculates the control inputs for the actual trajectory.

•
The MPC framework's RTG module is updated with the actual trajectory at the update point, which is before the IBP.

Validation of the Model Predictive Control Framework
The validation of the MPC framework is done by incorporating the controller in the Simulink ® model, using lvl-2 custom S-functions. A series of realistic telegraph trajectories, generated by a probability density function from historical data, was provided by Damen Naval BV for the vessel. Using the combinator curve presented in Section 3.1, these telegraph trajectories correspond to the disturbances for the propulsive load and the shaft speed. The hotel load power demand has been constant at 1000 kW.
The above loading sequences were spliced together to construct two different scenarios, referred to as Profiles 1-2. The profiles are used in an information barrier scheme which has been devised for the validation of the MPC framework: Two online optimal DP controllers are being fed the same disturbances with the MPC framework controller. At a certain time-step, called the information barrier point (IBP), the telegraph trajectory changes from the initial prediction to the actual trajectory (solid and dashed lines, respectively, shown for profiles 1 and 2 in Figure 7):

•
The first online DP controller has perfect information about the disturbances, including the change at the IBP. This means that it is the optimal controller and calculates the control inputs for the telegraph trajectory spliced at the barrier point.

•
The second online DP controller solves for the initial prediction until the IBP, where it recalculates the control inputs for the actual trajectory.

•
The MPC framework's RTG module is updated with the actual trajectory at the update point, which is before the IBP. At the final horizon, the MPC framework can either start diminishing its horizon or use the last calculated control inputs from its RTG module. The first option should be preferred. However, it has not been implemented here. For this reason, the final horizon area is not included in the online simulation, resulting in comparison problems when the optimal and barrier DP controllers have not used the same amount of battery charge at the end of the online simulation. Offline simulations with the reduced model descriptions are included to help with interpreting the results. First, the results for the CD mode will be presented, were the end SoC is free to terminate anywhere within [0. 2 1],and then the CS mode results will be presented, where the SoC starts at 0.9 and can only terminate with a value in [0.89 0.91]. Table 2 shows the results for CD battery mode. For profile 1, not all of the battery charge is depleted by the barrier controller, and for profile 2, no such effect takes place. The barrier benchmark controller optimally solves a control problem twice: the first time based on the full initially predicted trajectory up until the IBP at t = 2400 s, after which it solves a second time, now for the actual trajectory (remaining mission). The optimal controller solves only once for the actual trajectory. The MPC framework has an information advantage compared to the barrier controller due to the update point at t = 600 s, however, its optimality is not guaranteed and is being benchmarked. It is noted here that the consumption reduction of the MPC framework compared to the barrier controller is attributed to the information advantage described above. However, even if there is no such consumption reduction, the framework's performance should be further assessed for its robust behaviour (i.e., yielding a feasible control policy while satisfying all constraints) and for its consumption when compared to the optimal controller.
In Figure 8, in the plot for the online simulation of Profile 1, the MPC framework SoC trajectory (blue) starts to follow the optimal controller SoC trajectory (red) after the update point at t = 600s.This is the point where the RTG module of the MPC framework solves for the updated actual trajectory and provides the new reference to the MPC module of the framework. In the case of the MPC framework, the update for the trajec- At the final horizon, the MPC framework can either start diminishing its horizon or use the last calculated control inputs from its RTG module. The first option should be preferred. However, it has not been implemented here. For this reason, the final horizon area is not included in the online simulation, resulting in comparison problems when the optimal and barrier DP controllers have not used the same amount of battery charge at the end of the online simulation. Offline simulations with the reduced model descriptions are included to help with interpreting the results. First, the results for the CD mode will be presented, were the end SoC is free to terminate anywhere within [0. 2 1],and then the CS mode results will be presented, where the SoC starts at 0.9 and can only terminate with a value in [0.89 0.91]. Table 2 shows the results for CD battery mode. For profile 1, not all of the battery charge is depleted by the barrier controller, and for profile 2, no such effect takes place. The barrier benchmark controller optimally solves a control problem twice: the first time based on the full initially predicted trajectory up until the IBP at t = 2400 s, after which it solves a second time, now for the actual trajectory (remaining mission). The optimal controller solves only once for the actual trajectory. The MPC framework has an information advantage compared to the barrier controller due to the update point at t = 600 s, however, its optimality is not guaranteed and is being benchmarked. It is noted here that the consumption reduction of the MPC framework compared to the barrier controller is attributed to the information advantage described above. However, even if there is no such consumption reduction, the framework's performance should be further assessed for its robust behaviour (i.e., yielding a feasible control policy while satisfying all constraints) and for its consumption when compared to the optimal controller.
In Figure 8, in the plot for the online simulation of Profile 1, the MPC framework SoC trajectory (blue) starts to follow the optimal controller SoC trajectory (red) after the update point at t = 600 s.This is the point where the RTG module of the MPC framework solves for the updated actual trajectory and provides the new reference to the MPC module of the framework. In the case of the MPC framework, the update for the trajectory occurs so that the controller is able to deplete the battery in time, while the late informed barrier controller (yellow) finishes the mission with some "trapped" battery charge. The offline results in Table 2 and Figure 8 also show the effect clearly. tory occurs so that the controller is able to deplete the battery in time, while the late informed barrier controller (yellow) finishes the mission with some "trapped" battery charge. The offline results in Table 2 and Figure 8 also show the effect clearly. The effect of "trapped" battery charge can be understood in Figure 9, for Profile 1. At the region between 700 and 1000 s the barrier controller does not use the PTI, while the MPC framework, updated with the new information tries to use up the battery charge it did not use at the beginning of the mission. Not only that, but at 1500 s, the MPC framework depletes the battery more aggressively than the optimal controller to catch up. After 2400 s the propulsive load is too small for the PTI to be used, and the barrier controller cannot use any remaining charge.
Profile 2 presents little room for alternate control "paths", and a big portion of the optimal and barrier controllers coincide. The MPC framework follow very closely the setpoints of the optimal controller as it can be seen in Figure 9.
Two more profiles were tested for CD mode with similar results to Profile 2, with the MPC framework following consistently the optimal controller solution after the update point.
A parametric study for the position of the update point was carried out. Table 3 presents the results and demonstrates that if the "trapped" charge effect is present, there is a relatively wide update window for the MPC framework to yield savings and the savings are reducing consistently with later updates. In Figure 10 (right) it is shown that the PTO/I activity is consistent with the update points. For example, only the controller with an update point at 600 s follows the PTI setpoints of the optimal controller between 800 and 1000 s, while it can be seen that for the controller with an update point at 1800 s the highest PTI setpoints are observed after that 1800 s mark in an attempt to deplete the "trapped" charge. The effect of "trapped" battery charge can be understood in Figure 9, for Profile 1. At the region between 700 and 1000 s the barrier controller does not use the PTI, while the MPC framework, updated with the new information tries to use up the battery charge it did not use at the beginning of the mission. Not only that, but at 1500 s, the MPC framework depletes the battery more aggressively than the optimal controller to catch up. After 2400 s the propulsive load is too small for the PTI to be used, and the barrier controller cannot use any remaining charge.
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 20 of 28 A second parametric study compared different values for the hotel loads. It was found that there are fuel savings when the electric machines are used, either in PTI or PTO mode (Figure 10a-design hotel load 7269 kW). In the latter case the controller takes power off the shaft and offsets the main engine setpoints to the low consumption region (red cluster at 7MW on the main engine SFOC plot- Figure 11). For intermediate hotel loads (5 MW) the energy in the battery matches the hotel energy consumption of the mission and results in low PTO/PTI activity. The constraint on the torque of the electric machine is active and without the PTO/PTI mechanism to offset the setpoints of the main engine no reduction in fuel consumption was observed (Table 4).  Profile 2 presents little room for alternate control "paths", and a big portion of the optimal and barrier controllers coincide. The MPC framework follow very closely the setpoints of the optimal controller as it can be seen in Figure 9.
Two more profiles were tested for CD mode with similar results to Profile 2, with the MPC framework following consistently the optimal controller solution after the update point.
A parametric study for the position of the update point was carried out. Table 3 presents the results and demonstrates that if the "trapped" charge effect is present, there is a relatively wide update window for the MPC framework to yield savings and the savings are reducing consistently with later updates. In Figure 10 (right) it is shown that the PTO/I activity is consistent with the update points. For example, only the controller with an update point at 600 s follows the PTI setpoints of the optimal controller between 800 and 1000 s, while it can be seen that for the controller with an update point at 1800 s the highest PTI setpoints are observed after that 1800 s mark in an attempt to deplete the "trapped" charge. A second parametric study compared different values for the hotel loads. It was found that there are fuel savings when the electric machines are used, either in PTI or PTO mode (Figure 10a-design hotel load 7269 kW). In the latter case the controller takes power off the shaft and offsets the main engine setpoints to the low consumption region (red cluster at 7MW on the main engine SFOC plot- Figure 11). For intermediate hotel loads (5 MW) the energy in the battery matches the hotel energy consumption of the mission and results in low PTO/PTI activity. The constraint on the torque of the electric machine is active and without the PTO/PTI mechanism to offset the setpoints of the main engine no reduction in fuel consumption was observed (Table 4).  A second parametric study compared different values for the hotel loads. It was found that there are fuel savings when the electric machines are used, either in PTI or PTO mode (Figure 10a-design hotel load 7269 kW). In the latter case the controller takes power off the shaft and offsets the main engine setpoints to the low consumption region (red cluster at 7 MW on the main engine SFOC plot- Figure 11). For intermediate hotel loads (5 MW) the energy in the battery matches the hotel energy consumption of the mission and results in low PTO/PTI activity. The constraint on the torque of the electric machine is active and without the PTO/PTI mechanism to offset the setpoints of the main engine no reduction in fuel consumption was observed (Table 4). The profiles were also tested for CS mode. For the profiles in CS mode the utilisation of the battery was below 0.3% of the total installed charge and the induction machine was only initiated during low propulsive demand, where the operating limit of the main engines was activated. For Profile 1 the experiment was run again with the PTO/I torque limit relaxed from 0.3 to 0.1 and the hotel demand raised to 7269 kW. The results are shown in Figure 12, where it can be seen that the MPC framework closely follows the optimal controllers. It is noted that no fuel savings were observed in any of the CS modes, results consistent with the findings in [6] were CS fuel savings due to battery utilisation were not found.  The profiles were also tested for CS mode. For the profiles in CS mode the utilisation of the battery was below 0.3% of the total installed charge and the induction machine was only initiated during low propulsive demand, where the operating limit of the main engines was activated. For Profile 1 the experiment was run again with the PTO/I torque limit relaxed from 0.3 to 0.1 and the hotel demand raised to 7269 kW. The results are shown in Figure 12, where it can be seen that the MPC framework closely follows the optimal controllers. It is noted that no fuel savings were observed in any of the CS modes, results consistent with the findings in [6] were CS fuel savings due to battery utilisation were not found. The profiles were also tested for CS mode. For the profiles in CS mode the utilisation of the battery was below 0.3% of the total installed charge and the induction machine was only initiated during low propulsive demand, where the operating limit of the main engines was activated. For Profile 1 the experiment was run again with the PTO/I torque limit relaxed from 0.3 to 0.1 and the hotel demand raised to 7269 kW. The results are shown in Figure 12, where it can be seen that the MPC framework closely follows the optimal controllers. It is noted that no fuel savings were observed in any of the CS modes, results consistent with the findings in [6] were CS fuel savings due to battery utilisation were not found.

Discussion
The presence of a shipboard battery pack introduces loading sequence-specific, rather than profile-specific, system behavior. This should be addressed by control policies aiming for optimal fuel consumption and effective enforcement of the system constraints. The case study demonstrates the following about the proposed MPC framework:

•
The MPC framework shows close-to-optimal performance and satisfies all operating constraints, including the integral constraints on the SoC, for all tested loading sequences and modes. A time interval of 600 s for both the prediction and the control horizon has been found to be sufficient for achieving such performance.

•
In CD mode the consumption reduction due to utilisation of mission-scale predictions can be substantial. Compared to the fully-informed, optimal controller, the MPC framework selects higher power setpoints only after it receives the information for the updated profile, both in the standard and the parametric study (Figures 9a  and 10b). This means that the primary factor for reducing the fuel consumed is the use of all the energy in the battery. This in turn offsets the consumption of the main engines and the generator to lower values. A reduction of up to 3.5% was achieved for this specific case study in CD mode and under specific loading sequences. In CS mode, savings due to unused charge avoidance are not possible.

Discussion
The presence of a shipboard battery pack introduces loading sequence-specific, rather than profile-specific, system behavior. This should be addressed by control policies aiming for optimal fuel consumption and effective enforcement of the system constraints. The case study demonstrates the following about the proposed MPC framework:

•
The MPC framework shows close-to-optimal performance and satisfies all operating constraints, including the integral constraints on the SoC, for all tested loading sequences and modes. A time interval of 600 s for both the prediction and the control horizon has been found to be sufficient for achieving such performance.

•
In CD mode the consumption reduction due to utilisation of mission-scale predictions can be substantial. Compared to the fully-informed, optimal controller, the MPC framework selects higher power setpoints only after it receives the information for the updated profile, both in the standard and the parametric study (Figures 9a and 10b). This means that the primary factor for reducing the fuel consumed is the use of all the energy in the battery. This in turn offsets the consumption of the main engines and the generator to lower values. A reduction of up to 3.5% was achieved for this specific case study in CD mode and under specific loading sequences. In CS mode, savings due to unused charge avoidance are not possible.

•
When the induction motors are active, the reduction in fuel consumption due to mission scale information is possible (0.4% reduction observed). The control inputs from the MPC framework yield fuel-efficient operating points in the envelope of the main engine. This is achieved by storing energy, via the induction motors, to the battery plant and using it at a later time, in such a way that the fuel consumed is minimised.

•
For several loading sequences, the two benchmark solutions coincide partially or completely. In such cases, few solutions that satisfy all of the OCP's constraints exist, leading to a solution space where fuel consumption reduction due to use of mission scale information is insignificant. In the simulations for the MPC framework, the DP solvers in the RTG and MPC modules of the framework have been shown to be effective in finding these feasible solutions, adding to the controller's robustness.

•
The computational load allows for an embedded controller. The controller framework implemented in Simulink ® , simulated together with the dynamic model, require less time to compute than the mission on a laptop computer.
The results of the study show the effectiveness of the proposed framework, compared to the state-of-the-art shipboard EMSs, particularly when handling mission-scale effects and integral constraints introduced by the battery system. Furthermore, the performance of the proposed framework in terms of fuel consumption and robustnessis close to the optimal solution. More specifically, it is demonstrated that the controller yields feasible control policies and is robust regarding the integral constraint on the battery SoC, addressing concerns identified in the literature, such as the insufficient peak demand availability observed in [25]. It is also capable of adapting to real-time updates on the future disturbance estimations, as opposed to the controller proposed in [38]. At the same time, it addresses the problem of instantaneous EMS (ECMS, A-ECMS), identified in [16], where the performance of the controller is poor under off-design power demand sequences, different from those used to tune the controller's equivalency factor(s).
These results may motivate future work to focus on incorporating online short-term and long-term prediction schemes into the MPC framework, in an attempt to design a complete controller solution (see Section 1.3). The lack of operating profile-specific parameters may enable early ship power plant design studies to quantify fuel consumption savings and battery sizing for different operating profiles and sequences. Finally, complex ship power plants can benefit from the framework's ability to readily incorporate additional control input variables and states as well as terms in the OCP's objective function. Acknowledgments: Special thanks goes to Damen Naval BV for supplying the Simulink ® model of the naval vessel and providing necessary information and support throughout this study.

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

Appendix A
In order to find how to divide the total power demand among the available generator sets, the active number of generator sets can be found by an integer division NoE = div(P gen ,P genset,nom ) + 1, and then, by dividing the total power with this number the consumption can be found: The above yields optimal results in the common case where (a) the fuel consumption of the generator set is a convex function of its power demand and (b) the generators are identical.
By considering two generator sets with power setpoints P g,1 ,P g,2 , what needs to be proven is that if point C is the common operation point of the two engines, there exist no points A and B, that cover the total power demand, such as: m f ,c ⇔ f (P a,g,1 ) + f (P b,g,2 ) < 2 f (P c,g,12 ), P a,g,1 < P c,g,12 < P b,g,2 This can be proven from the definition of convexity: Now, since point C is between A and B, it can be expressed in the form: P c,g,12 = tP a,g,1 + (1 − t)P b,g,2 (A4) which yields: t = P c,g,12 − P b,g,2 P a,g,1 − P b,g,2 (A5) and by substituting in the equation above the expression 2P c,g,12 = P a,g,1 + P b,g,2 = P gen , which means that in both scenarios, the total power has to add up to the total power demand P gen , yields that t = 1/2. From the definition inequality (A3), and by substituting the m f function, it is now: which can be generalised for more engines and which proves that the common power setpoint P c,g,12 = P gen,i from Equation (A1) is optimal. In the case of different engines, but their consumption functions are still convex, the formulation above can be extended, however, that case is not discussed here. Regarding the operating limits, the following inequalities should hold: P gen,i ∈ [P gen,min , P gen,max ] P gen ≤ k g P gen,max where k g the number of generator sets available on the vessel.