Open access peer-reviewed chapter

An Adaptive Load Frequency Control Based on Least Square Method

Written By

Adelhard Beni Rehiara, Naoto Yorino, Yutaka Sasaki and Yoshifumi Zoka

Submitted: April 11th, 2019 Reviewed: October 28th, 2019 Published: April 1st, 2020

DOI: 10.5772/intechopen.90300

Chapter metrics overview

746 Chapter Downloads

View Full Metrics


Modern power system becomes a complex system consisted with various load and power stations. Therefore, it may spread into some areas of power system in neighborhood, and so a load frequency control (LFC) is a necessity device to regulate the frequency of the power system by distributing the load to the generating units and controlling tie-line power interchange between areas. The integration of renewable energy sources (RES) into a power grid has presented important issues about stability and security of power system. In such conditions, the use of conventional LFC may not be sufficient to protect the power system against the power changes. In this chapter, an adaptive LFC controller is proposed using the least square method (LSM). The controller adopts an internal model control (IMC) structure in two scenarios, i.e., static controller gain with adaptive internal model and both the adaptive controller gain and adaptive internal model. A two-area power system is used to test and to validate both performance and the effectiveness of this controller through some case studies.


  • adaptive controller
  • power system stability
  • LSM
  • LFC
  • IMC
  • res

1. Introduction

Demand load in a power system is continuously varied by the time and the change of the active and reactive power demand has introduced generator-load mismatching power. When the load increases, it will slow down the rotor speed which result drop of the frequency. On the other hand, when the load is reduced, the load frequency will rise. The change of load frequency will directly affect electrical motor performance, further interfering the protection of the system [1].

Smart grid technology is well developed, and it currently has been widely used in the power system operation due to the need for integrating renewable energy to the existing power grid. Penetration of the renewable energy resources such as solar generations and wind turbines to the power grid has introduced significant issues of stability and security of the power system. However, frequency stability is a major issue in the power system operation due to continuous output change of the renewables. Therefore a load frequency control (LFC) is an essential device to back up the automatic generation control (AGC) to keep the frequency stable during power system operation [2].

The frequency controllers in a power system consist of AGC and LFC systems as primary control and secondary control. An AGC will respond to small changes of frequency in a generator, and an LFC will have to regulate load frequency in a large area of the power system and large change in load frequency which is typically in between ±5 and 6% of the frequency bias [2, 3, 4, 5].

In this chapter, adaptive LFC controllers are introduced, where an adaptive IMC model that is repeatedly updated by the least square method (LSM) in real-time operation is proposed. It is shown that the target model is successfully identified and therefore that the proposed LFC controller scheme effectively keeps the system frequency at a desired set point. The effectiveness of the proposed controller is demonstrated by the simulations using a standard LFC model representing two-area interconnected power system.


2. Power system model

2.1 Mathematical modelling of generator

Eq. (1) is well known as a swing equation that describes the rotor dynamics and hence is known as the swing Equation [6, 7]. The internal EMF angle δ is called the load angle that indicates how much power can be transferred:


For small perturbation, the swing equation of a synchronous machine will be formed in Eqs. (2) and (3) for its small deviation in speed. Therefore, Laplace transformation of Eq. (3) is shown in Eqs. (4) [6, 7]:


The mathematical model of the generator as formulated in Eq. (4) can be figured in Figure 1.

Figure 1.

Block diagram for generator.

2.2 Mathematical modelling of load

The load on a power system consists of a variety of electrical devices which is a resistive or inductive load. The equipment used for lighting purposes and heating are basically resistive in nature, and this kind of load is independent to frequency. On the other hand, rotating devices, such as fans and pumps, are basically a composite of the resistive and inductive components which are dependent to the frequency changes. The speed-load characteristic of the composite load is given by [6, 7]:


where ΔPL is the load change and DΔf is the frequency sensitive load change. D is expressed as percent change in load by percent change in frequency. By adding the load to the generator, block diagram of both load and generator is figured out in Figure 2.

Figure 2.

Generator and load block diagram.

2.3 Mathematical modelling for prime mover

The electrical energy is generated inside a power generation by converting the other kind of energy sources by means of a prime mover. The prime mover may be diesel machines, hydraulic turbines at waterfalls, or steam turbines. The model for the turbine relates the changes in mechanical power output ΔPm to the changes in the steam valve position ΔPV [6, 7]:


Figure 3 expresses the prime mover block diagram, where Tt is the turbine constant which has the range in between 0.2 and 2.0 seconds.

Figure 3.

Block diagram for prime mover.

2.4 Mathematical modelling for governor

The electrical power will exceed the mechanical power input when the electrical load is suddenly increased. This condition will result in the extraction from the rotating energy of the turbine. Then the kinetic energy stored in the machine is reduced and slows down the speed of prime mover. Therefore, to compensate the reduced speed, the governor sends a command to supply more volumes of water or steam or gas to increase the prime mover speed.

Speed regulation R is given as the curve slope in Figure 4. The typical speed regulation values of generator are in between 5 and 6% from zero to maximum of load [6, 7]:

Figure 4.

Speed drop regulation.


The quantity of ΔPg is converted to the position of steam valve ΔPV by a governor time constant Tg. Therefore, the s-domain relation of ΔPV and ΔPg is a linear relationship by considering the simple time constant Tg [6, 7] (Figure 5):

Figure 5.

Block diagram for governor.


Finally, Figure 6 summarizes a combining of all of the block diagrams from earlier block diagrams for a single area system.

Figure 6.

Completed power system block diagram.

A completed LFC block diagram for multi-area power system, including controller, frequency bias, and tie-line power change, can be redrawn in Figure 7 [2, 3, 4, 8].

Figure 7.

Power system dynamics.

The tie-line power change Ptie is calculated for all area n using Eq. (9), and the area control error (ACE) which is a suitable linear combination of frequency f and tie-line power changes for each area is found using Eq. (10) as follows [3, 4, 9]:


A general state-space model is used to describe the power system model as shown in Eqs. (11) and (12):


where x(t), u(t), w(t), and y(t) are the matrices of state variables, input variables, control variable, and output variable, respectively. Four variables of the state variables are ΔPg,i, ΔPm,i, Δfi, and ΔPtie,I, and the input variables are ΔPL,i and Δvi. ΔPc,I is the control variable, while ACEi is the output variable.

Due to no direct connection between input and output variables, the feed forward matrix D is removed from the model. Therefore the system matrices of a LFC system are written in Eqs. (13)(16) [3, 4, 9]:


where Pg,i, Pm,i, PL,i, Pc,i, yi, Hi, di, Ri, βi. Tij, Tg,i, and Tt,i are the output of governor, the prime mover power, the load, the control action, the output of system, the inertia constant, the damping coefficient, the characteristic of speed droop, the bias factor of frequency, the tie-line synchronizing coefficient between reference area i and area j, and the time constant of governor and turbine, respectively.

2.5 System response of power change

Consider a single machine system connected to an infinite bus as shown in Figure 8, and its swing equation in steady-state condition can be expressed in Eq. (17) [6]:

Figure 8.

Single machine system.


If there is some change in mechanical power input Pm as the result of disturbances or load changes, power angle δ will change to a new state as δ = δ0 + Δδ. Then it will further influence swing equation into Eq. (18). The change affects the swing equation in terms of incremental changes in power angle as in Eq. (19):


The quantity of Pmax cos δ0 is known as the synchronizing coefficient Ps which is the slope of power angle curve at δ0. The root(s) of the second-order differential equation in Eq. (19) can be shown in Eq. (20):


It can be seen from Eq. (20) that there are possibilities of roots in s-plane when Ps is either positive or negative. A root in the right hand s-plane, where causes system unstable and responses increased exponentially, is gotten when synchronizing coefficient Ps is negative. In other way, two roots will be on j-ω axes of s-plane for Ps negative that causes system responses, oscillatory and undamped with natural frequency as in Eq. (21):


3. Controller structure

3.1 Model predictive control

Model predictive control (MPC) is an advance optimal control in the field of control systems engineering. In an MPC, the optimal trajectory movement is given by properly choosing the MPC gain so that the control errors can be minimized.

The objective of the predictive control is to compute the manipulated variable u in order to optimize the output behavior of a controlled plant y. An MPC will use its internal model to calculate the manipulated variable [10] (Figure 9).

Figure 9.

Principle of MPC.

At a given discrete time k, the plant output is estimated through prediction horizon Np from time k + 1 to k + Np, and the MPC controller output is predicted by control horizon Nc. The output of the plant is continued to be minimized based on specified objective function which is typically in the form of a quadratic function as shown in Eq. (22) [10]:


where N, Q, and R are the error weight matrices of the control action, output system, and rate of change in control action, respectively. Δŷ(k + m|k) is the output of plant in prediction, ∆u(k + m|k) is the change rate in control action under prediction, and u(k + m|k) is the prediction of the action of optimal control, for the measurements given at time k + m from the reference time k.

In order to have optimal result, prediction horizon and control horizon have to be set properly so that the MPC controller can work in high performance. On the other hand, absence of doing this will cause the MPC to lose optimal action that will result in a high overshoot response.

In a short time control horizon, Nc will respond with high control action that results in an overshoot after the end of control horizon time. While long time control horizon bring the controller be aggressive and used much energy to accelerate and decelerate the control action in order to keep it constant at set point. In same way, too long prediction horizon Np will drop controller performance due to the extra time needed for calculating the trajectory movement, while short time prediction horizon will cause inaccurate calculation of the trajectory.

For a given a linear system in continues time, a system state condition is expressed in Eq. (23) while Eq. (24) shows its formulation for discrete system with sampling time k.


Then the model will be converted to an augmented model so that later the quadratic programming problem with respect to ∆U could be formed easily. Because u(k) = u(k-1) + Δu(k), then Eq. (24) can be rewritten as in Eq. (26), and its state space from is given in Eqs. (27) and (28) [11]:


For the Np prediction horizon based on the above equation model, the output of the prediction could be written as:


Eq. (29) can be simplified as in Eq. (30):


Then the objective function could be written as:


Therefore, Eq. (31) could be solved efficiently by a quadratic programming.

3.2 Internal model control (IMC)

An IMC controller structure can be seen on Figure 10 where an internal model is used parallel with the plant. This internal model, which is also called efferent model, used plant model to estimate plant output. It is known that to find the internal model as same as plant model is difficult due to plant dynamics that is not well captured in the modeling part.

Figure 10.

IMC structure.

In an IMC controller, an efferent model is expected to correct the actual output before it is feedback to the controller. This is the difference between MPC controllers with classical controllers. In order to get maximum performance, the efferent model may have input and output relation as close to the real plant model so that the only feedback signal is the disturbance.

As a type of controller, an IMC can control a plant directly. This controller can also be used to tune another controller if it is difficult to be the best parameter for this controller. Based on the IMC structure, this main controller can be any type of controller so that this IMC can be combined with any type of controllers [12]. The control law for an IMC control is written in Eq. (32) to Eq. (34) [13, 14]:


The difference between classical controller and an IMC is that an IMC will correct the actual output before it is fed back. Since an IMC uses the efferent model, the model should be a perfect model to have the highest control performance. The way to provide the model in an IMC can be in forward model, inverse models, combination of both forward and inverse models, or adaptive model.

3.3 Adaptive model structure

Due to highly nonlinearity in a power system, the classical model of power system may not be accurate to configure the real power system. Therefore the other objective of this paper is to find simple model of a load frequency of power system using the least square method (LSM). In order to capture the essential dynamics of power system, a first-order lag model is adopted in this chapter as follows:


The power system model, including governor, turbine, and rotating mass and load model as in Figure 7 will be replaced by the first-order lag model in Eq. (35). Therefore it has become a simple model as formulated in Eq. (36):


Expanding Eq. (35) gives


For discrete time system, the above equation can be described in a matrix form as follows:


Consider a linear equation of least square method y = Hx; a solution for minimizing the error can be written as [15, 16]:


Then expanding J(x) gives


Taking the derivative for the J(x) gives


Minimizing the derivative by setting it to zero gives


If HTH is invertible, then the least square solution is given in Eqs. (44) and (45) for its simple form:


The least square solution for the discrete time system of the power system model is shown in Eq. (46).


where H is the pseudo inverse of the H matrix. A and B matrices are taken from the LFC state-space model, and ACE is selected as input u to the LFC controller.

3.4 Model simplification

In order to adopt the least square solution into the internal model of the LFC system, the complex system of the LFC is likely to be reduced into a first-order model. After the disturbance is entering the LFC system, generations are expected to change as fast as possible to meet the load demand. This expectation may not be achieved due to slow response of governor and turbine operations. For these purposes, both turbine and generator responses are neglected to derive a simple expression for the time response.

Considering an LFC model in Figure 11, the model in steady-state condition at ΔPL = 0 will be a third-order transfer function as written in Eq. (47):

Figure 11.

Power system model.


where Kg, Kt, and Kp are a constant of the governor, turbine, and power system, respectively.

Typically in a LFC system, power system time constant is relatively higher than governor and turbine time constant as described in Eq. (48). Therefore both governor and turbine time constant are negligible, and by adjusting KgKt = 1, Eq. (47) is simplified to be a first-order system in Eq. (49):


The first-order system effectively reduced the unexpected response. An example of system responses of first-order and third-order systems is given in Figure 12.

Figure 12.

First- and third-order system responses.

The dynamics response of a frequency control means how the frequency is immediately corrected after the disturbance and before the system reaches new steady-state condition [17]. In this case, a static first-order model may not solve the LFC problem since the disturbance is fluctuated by the time. Therefore an adaptive model of first-order system is suggested to be used in solving this problem.

3.5 Computational procedure

Power system model is expressed in a first-order lag system. The way to build the model can be a combination of either inertia constant and damping or a gain and time constant. In this case, the second combination will be used to provide the power system model using the least square method. The power system equation can be written as follows.

Block diagram of the adaptive LFC system is shown in Figure 13 as rewritten from [1]. The adaptive internal model is utilized to provide an update model of the power system, and the model is generated by using a LSM method. The proposed adaptive LFC controller uses MPC controller as its main controller combining with an adaptive efferent model as the internal model.

Figure 13.

Adaptive LFC system using IMC structure.

In order to reach the best controller solution, the simulation will be done using the given algorithms. IMC type 1 includes the adaptive internal model, while in IMC type 2, controller gain Kmpc is also updated with the change of internal model. Therefore, IMC type 2 uses an algorithm of both adaptive internal model and adaptive MPC controller.

Algorithm 1:
IMC with adaptive model (type 1)
Algorithm 2:
IMC with adaptive model and controller (type 2)
1.   set disturbances and noises
2.   configure presetting model
3.          ⁝
4.   for j = 1 to simulation time
5.         ⁝
6.       for i = 1 to n-area
7.            calculate ∆u
8.            upgrade state matrix ∆
9.           if j = desired time then
10.                 update K & T by LSM
11.                 change internal model
12.           end
13.           calculate model output ∆ym
14.          correct signal ∆ = ∆y − ∆ym
15.           ⁝
16.      end
17.  end.
1. set disturbances and noises
2. configure presetting model
3.       ⁝
4. for j = 1 to simulation time
5.      ⁝
6.     for i = 1 to n-area
7.          calculate ∆u
8.          upgrade state matrix ∆
9.           if j = desired time then
10.             update K & T by LSM
11.             change internal model
12.            define new MPC gain Kmpc
13.         end
14.         calculate model output ∆ym
15.          correct signal ∆ = ∆y − ∆ym
16.          ⁝
17.    end
18. end.


4. Controller test

The investigated power system consists of a two-area power system which is modified from a three-area poser system in [2, 3, 4, 9, 18]. The controller performance is tested by comparing a classical MPC controller and the proposed adaptive IMC controller. System configuration is based on Figure 7 where the system parameter is shown in Table 1.

The proposed adaptive LFC controller type 1 and type 2 will be computed as follows. The gain K from PC to Δf in the detailed model is set to 66.5, which is presetting the value for each area. An initial value of the time constant T is chosen as the same value as the power system model of Figure 14, which is set equal to 1. Those will be updated online based on system identification method starting from t = 20.

Figure 14.

Multi-area power system configuration.

The main controller for the proposed controller is an MPC controller which at initial condition is set as the same as the existing controller. Prediction horizon Np = 10 and control horizon Nc = 2 are applied to both existing and proposed controllers. In the whole of simulation period, the existing MPC controller gain is not change where computed off line in the beginning before simulation is started.

Three different cases are performed in simulations under disturbance setting as in Table 2. The disturbances are load changes which are imposed by white noises. The random disturbance is white noises with a maximum magnitude about 0.05 pu, while the step disturbance is a sudden change of load. The step disturbance is assumed as the action of economic load dispatching which is applied at t = 40 for case II and case III in addition to the white noise.

4.1 Simulations

A. Case 1.

Simulation results are shown in Figure 15 for Δf responses and in Figure 16 for ΔPm responses. In this case, only random disturbance is applied. It is observed that the adaptive controllers both type 1 and type 2 show slightly better performances compared to the existing MPC:

Figure 15.

Case 1: Δf responses.

Figure 16.

Case 1: ΔPm responses.

B. Case 2.

The step disturbance about 0.2 pu is applied at t = 40 s in addition to the random disturbance. The results are shown in Figure 17 for Δf responses and in Figure 18 for ΔPm responses.

Figure 17.

Case 2: Δf responses.

Figure 18.

Case 2: ΔPm responses.

A better performance is observed. It is noted that although the adaptive controllers are a much simpler configuration compared with the existing controller, the control performance is even better. This is an advantageous feature of the proposed method.

C. Case 3.

In this case, step disturbance about −0.1 pu. which is applied at t = 40 s assumed that load is released at t = 40 s including the random disturbance. The results are shown in Figure 19 for Δf responses and in Figure 20 for ΔPm responses.

Figure 19.

Case 3: Δf responses.

Figure 20.

Case 3: ΔPm responses.

The similar performance is shown. It is noted that the adaptive controllers can handle incoming disturbance and also released disturbance with even better performance compared to the existing controller.

4.2 Evaluations

Figures 2124 show how the internal model parameters are identified. It is observed from Figures 21 and 23 that initial gain K is updated as soon as the model identification process is completed at t = 20, which are consistent values based on the LSM. In the same way, Figures 22 and 24 show that the initial time constants T are updated very slightly around 1.0. Those values are continuously updated around the converged values.

Figure 21.

Gain K of IMC type 1.

Figure 22.

Time constant T of IMC type 1.

Figure 23.

Gain K of IMC type 2.

Figure 24.

Time constant T of IMC type 2.

LFC system based of IMC-1 is equivalent to the existing MPC controller in some cases, while the performance of IMC type 2 is even better in all areas and all cases. Tables 3 and 4 list measured values of the overshoot for the step disturbance and the standard deviations of both frequency and prime mover responses for all cases. It is seen from the tables that the performance of the proposed adaptive LFC controllers is equivalent or better than the conventional MPC controller. This implies that the adaptive LFC controllers can successfully identify the target model and handle the power system disturbances. In the same way, the controllers keep the system conditions successfully at the set points.

[pu s]
10.0160.20172.730.060.440.3827T12 = 0.20
20.0150.12472.820.070.300.3692T21 = 0.20

Table 1.

System parameters.

RandomStep change
1AppliedNot applied
2Applied0.2 pu.
3Applied−0.1 pu.

Table 2.

Disturbance setting.

Case 1Case 2Case 3
Area 1Area 2Area 1Area 2Area 1Area 2
Standard deviationMPC0.01860.01190.00940.00620.00330.0022

Table 3.

Frequency responses.

Case 1Case 2Case 3
Area 1Area 2Area 1Area 2Area 1Area 2
Standard deviationMPC0.09450.00340.04660.00180.00180.0006

Table 4.

Prime mover responses.

The simulations are carried out on PC with Intel Core i7 1.8 GHz CPU and 4 GB RAM using MATLAB 2016a under Windows 10. CPU times for the computation of controllers are listed in Table 5. It is observed that the adaptive LFC controllers are slower than the conventional method.


Table 5.

Simulation time (s).

In today’s power system condition, where smart grid and energy management system are applied, the power system dynamics is varying. Therefore an adaptive model would be an advantage and is required for the better performance of LFC system.

It is observed that time consuming is increased together with the improvement of the controller performances. Therefore the IMC type 2 controller needs more time to compute and update both internal model and controller gain while giving a better performance. On the other hand this time consuming is the order of millisecond and so it is acceptable for an LFC system which in general operated in the order of second. Therefore the proposed controller has bit complexities in hardware for computing the model and updating the controller since its consuming time is about 40% higher than the existing MPC controller.


5. Conclusions

This chapter presents adaptive LFC methods based on IMC controller structure, where the internal model is adaptively updated online in IMC type 1, while both internal model and MPC controller gain are restructured in IMC type 2 by using the least square method. The performance of the controller is fair in handling load disturbances in spite of relatively slow control cycle with ramp rate constraints of actual systems.

Simulation results show that the gain and time constant of the internal model have been adaptively changed. This change has guaranteed the better performance of the proposed controller. Based on the system responses, the adaptive IMC controller type 1 has responses similar to the existing MPC controller, while the adaptive IMC controller type 2 has shown its superiority compared to MPC controller and IMC type 1 controller. In contrary, consuming time becomes larger by the enhancement of the controller performance.


  1. 1. Rehiara AB, Chongkai H, Sasaki Y, Yorino N, Zoka Y. An adaptive internal model for load frequency control using extreme learning machine. TELKOMNIKA. 2018;16:1-6
  2. 2. He C, Rehiara AB, Sasaki Y, Yorino N, Zoka Y. The Application of Laguerre Functions Based Model Predictive Control on Load Frequency Control. In: the International Conference on Electrical Engineering (ICEE2018). Seoul, Republic of Korea; 2018. pp. 1190-1195
  3. 3. Rehiara AB, Sasaki Y, Yorino N, Zoka Y. A performance evaluation of load frequency controller using discrete model predictive controller. In: 2016 International Seminar on Intelligent Technology and Its Applications. Lombok, Indonesia; 2016. pp. 659-664
  4. 4. Rehiara AB, Chongkai H, Sasaki Y, Yorino N, Zoka Y. An adaptive IMC-MPC controller for improving LFC performance. In: 2017 IEEE Innovative Smart Grid Technologies—Asia. Auckland, New Zealand; 2018. pp. 1-6
  5. 5. Abdillah M, Setiadi H, Rehiara AB, Mahmoud K, Farid IW, Soeprijanto A. Optimal selection of LQR parameter using AIS for LFC in a multi-area power system. Journal of Mechatronics, Electrical Power, and Vehicular Technology. 2016;7(2):93
  6. 6. Kundur P. Power System Stability and Control. USA: McGraw Hill; 1993
  7. 7. El-Hawary ME. Electrical Power Systems. New Jersey: Prentice-Hall; 1983
  8. 8. Mohamed TH, Morel J, Bevrani H, Hiyama T. Model predictive based load frequency control-design concerning wind turbines. International Journal of Electrical Power & Energy Systems. 2012;43(1):859-867
  9. 9. Mohamed TH, Hassan AA, Bevrani H, Hiyama T. Model predictive based load frequency control design. In: 16th International Conference on Electrical Engineering. 2010
  10. 10. Wong C, Mears L, Ziegert J. Dead time compensation for a novel positioning system via predictive controls and virtual intermittent setpoints. In: 2009 International Manufacturing Science and Engineering Conference. Indiana, USA; 2009. pp. 1-8
  11. 11. Wang L. Model Predictive Control System Design and Implementation Using MATLAB. London: Springer-Verlag; 2009
  12. 12. Rehiara AB, Yorino N, Sasaki Y, Zoka Y. A novel adaptive LFC based on MPC method. IEEJ Transactions on Electrical and Electronic Engineering. 2019;14(8):1145-1152
  13. 13. Shigemasa T, Yukitomo M, Kuwata R. A model-driven PID control system and its case studies. Proceedings of International Conference on Control Applications. 2002;1:571-576
  14. 14. Jin Q, Feng C, Liu M. Fuzzy IMC for unstable systems with time delay. In: 2008 IEEE Pacific-Asia Workshop on Computational Intelligence and Industrial Application. Washington DC, USA; 2008. pp. 772-778
  15. 15. Skala V. Least square method robustness of computations: What is not usually considered and taught. In: Proceedings of the 2017 Federated Conference on Computer Science and Information Systems, FedCSIS. Prague, Czech Republic; 2017. pp. 537-541
  16. 16. More AN, Kohli PS, Kulkarni KH. Simple linear regression with least square estimation: An overview. International Journal of Computer Science and Information Technologies. 2016;7(6):2394-2396
  17. 17. Sivanagaraju S, Sreenivasan G. Power System Operation and Control. 1st ed. Delhi: Pearson Education; 2010
  18. 18. He C, Rehiara AB, Sasaki Y, Yorino N, Zoka Y. Model predictive load frequency control using unscented Kalman filter. In: IEEJ Technical Meeting on Power Engineering. Jeju Island, Republic of Korea; 2018

Written By

Adelhard Beni Rehiara, Naoto Yorino, Yutaka Sasaki and Yoshifumi Zoka

Submitted: April 11th, 2019 Reviewed: October 28th, 2019 Published: April 1st, 2020