Saïd Mammar

Naïma Ait Oufroukh

Zedjiga Yacine

Dalil Ichalal

Lydie Nouvelière

Invariant Set Based Variable Headway Time Vehicle Longitudinal Control Assistance

This paper presents the design and simulation test of an assistance system that helps the driver to perform vehicle following maneuvers. When activated, the assistance system drives the vehicle to the desired inter-vehicle distance. Based on Lyapunov theory and invariant sets, the assistance system is able to guarantee that the state trajectories remain bounded under bounded disturbance inputs, such as leading vehicle accelerations. The assistance system is designed in such a way that variable time headway strategy is naturally handled. The control law is computed using Linear and Bilinear Matrix Inequalities (LMI-BMI) methods. Some design parameters can be adjusted to handle the tradeoff between safety constraints and comfort specifications.

I. INTRODUCTION

Vehicle longitudinal control has been the object of intensive research over the last decades [START_REF] Lingyun | A comprehensive review of the development of adaptive cruise control systems[END_REF]. Starting with cruise control, which is now widely implemented in today's commercial vehicles, research has shifted to the adaptive cruise control (ACC) which is an extension of the cruise control that manages both relative speed and distance to a preceding vehicle [START_REF] Milanes | Comparing Fuzzy and Intelligent PI Controllers in Stop-and-Go Manoeuvres[END_REF]. Starting just before the new millennium, ACC is already available in luxury passenger cars and trucks. Several technological requirements have been also investigated in the context of automated highways, including sensor range characteristics, actuator performance and cooperation spacing policies. ACC or longitudinal control means acting on the throttle as well as the brake system. As mentioned in [START_REF] Naus | Design and implementation of parameterized adaptive cruise control: An explicit model predictive control approach[END_REF], ACC systems typically consist of two parts: a vehicle-independent part and a vehicle-dependent part [START_REF] Moon | Design, tuning and evaluation of a fullrange adaptive cruise control system with collision avoidance[END_REF], [START_REF] Prestl | The BMW active cruise control ACC[END_REF]. The vehicle dependent part in the inner loops directly controls the actuator, in order to achieve the desired acceleration/deceleration generated by the vehicleindependent part. Concerning this latter aspect, several strategies have been studied. In [START_REF] Lingyun | A comprehensive review of the development of adaptive cruise control systems[END_REF], different spacing policies have been cited in order to enhance both safety and capacity. In [START_REF] Jairo-Martinez | A Safe Longitudinal Control for Adaptive Cruise Control and Stop-and-Go Scenarios[END_REF], vehicle longitudinal control has been considered using a reference model inspired from a non-linear damper/spring model associated with safety zones. In [START_REF] Nouvelière | Experimental vehicle longitudinal control using a second order sliding modes technique[END_REF], five assistance modes for shared driving have been established and classified from the least constraining (instrumented mode) to the most restricting (automated mode) against the driver.

Time headway based spacing appears as the most promising strategy. It is actually implemented on commercial vehicles equipped with ACC, and it has been proved that it Authors are with Université d' Évry Val d'Essonne, France. IBISC: Informatique, Biologie Intégrative et Systèmes Complexes -EA 4526, 40 rue du Pelvoux CE1455, 91020, Evry, Cedex, France, (e-mail: (said.mammar, naima.aitoufroukh, zedjiga.yacine, dalil.ichalal, lydie.nouveliere)@ibisc.univ-evry.fr).

is natural for the driver. In fact a time headway of about 2 seconds has been added as a vehicle following requirement in driving rules. However, in steady state, time headway is directly related to lane capacity. Thus, from the road manager point of view, if the control of the vehicle is delegated to the ACC, it would be appreciable that time headway could be tuned in order to achieve macroscopic traffic characteristics requirements.

This paper focuses on the vehicle-independent part. It aims at developing a controller for a follower vehicle that achieves control requirements of comfort and safety for a given range of headway times. When writing down the model, the leading vehicle acceleration enters as an unknown bounded disturbance input. The variable is hard to estimate from available car industry sensors. The proposed assistance is aimed to ensure both vehicle following and collision avoidance in case of driver inattention. For this purpose, invariant sets framework, which has significantly developed in control engineering over the last decades, is used [START_REF] Blanchini | Set-theoretic methods in control[END_REF]. The existence of an invariant set (ellipsoidal or polyhedral set), under bounded disturbance, is equivalent to the "inputto state stability" (ISS) notion [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF]. This framework ensures that even if transitions exist between driver and assistance control of the vehicle, every trajectory that starts inside the invariant set will not exceed it, i.e. the resulting hybrid system trajectories will be bounded inside this set [START_REF] Minoiu | Driver steering assistance for lane departure avoidance based on hybrid automata and on composite Lyapunov function[END_REF].

The remainder of this paper is divided as follows. The assistance specifications are given in Section II including the description of the safety driving zones and a switching strategy between the driver and the assistance. Section III contains the description of the vehicle inter-distance model. The design of the longitudinal control assistance is further conducted. Section IV contains simulation results for several maneuvers including ACC and stop-and-go traffic. Conclusions are presented in Section V.

II. LONGITUDINAL CONTROL CONCEPT

A. Control objectives

The main goal of the proposed assistance is to achieve a kind of friendly and safe vehicle following. More specifically, the assistance should ensure vehicle following according to a chosen inter-distance policy and it should avoid rearend collisions due to a driver's lapse of attention during "normal following". The "normal following" means a driving situation during which the driver is following a leader vehicle without performing any special maneuver (e.g. hard braking, lane change), as well as keeping a certain distance to the leading vehicle. This will be defined in more details in the next section.

To achieve the above control goal, a vehicle following control law has been designed such that: (a) The closed loop system, corresponding to the vehicle following problem has to be asymptotically stable to zero steady state. The reference distance could be based on variable time headway in order to allow the lower or greater lane capacity when required by the infrastructure manager. (b) During the assistance intervention, the vehicle should not enter the collision zone. (c) The state variables and the vehicle acceleration/brake have to be bounded to guarantee safety and comfort.

B. Safety zones related to the car following

In order to characterize different safety levels, the following zones have been defined (see Fig. 1(a) and 1(b)):

• In the normal green zone, the vehicle is in a safe operating region. The spacing distance s is around the nominal inter-distance. The assistance should switch on if the driver is not attentive. • The risky orange zone, incorporates two ranges. In the first one (see region [START_REF] Bando | Analysis of Optimal Velocity Model with Explicit Delay[END_REF] in Fig. 1(a)), the spacing is lower than the reference spacing. The assistance should override the driver in order to avoid reaching a collision red zone and to help him return in the green zone.

Reciprocally, in the second orange range (see region [START_REF] Blanchini | Set-theoretic methods in control[END_REF] in Fig. 1(a)), the inter-distance is higher than the reference distance. The assistance could also override the driver and accelerate in order to maintain the car following mode, if the requested acceleration is feasible.

If not, a cruise control mode could be activated. • In the collision red zone, the spacing is under the bias safe distance and a collision might occur if the leading vehicle performs a hard braking. In the other case, the cruising red zone, the inter-distance is too high and the vehicle is no more in a car following mode. The assistance deactivates and a cruise control mode takes its place. All the variables presented in Fig. 1(a) and 1(b) are described in Sections III-A, III-F and III-G.

C. Switching strategy

The switching strategy for assistance activation or deactivation uses of an invariant set ε p which includes the "normal following" zone and is located inside the limits defining the "safe following" zone (see Fig. 1(a)). Switching occurring inside this invariant set ensures that the controlled vehicle dynamics are stable and remains within this invariant set. Activation regions are shown in Fig. 1(a) and the switching follows the rules described below:

• The assistance is always activable inside the "normal following" region when the driver is not attentive. Outside the "normal following" region, the switching on occurs in the following context:

• When the leading vehicle is too far and approaching, the assistance is activable, but only inside the invariant set ε P . • When the leading vehicle is too far and going away, the assistance is not activable. The vehicle should enter in the cruise control mode. • In case the leading vehicle is close and going away, the assistance overrides the driver only if the system state is in the invariant set ε P . • When leading vehicle is too close and approaching, the assistance overrides the driver while the operating point is still inside the invariant set ε P .

III. VEHICLE MODEL AND CONTROL SYNTHESIS

A. Vehicle model

Vehicle longitudinal dynamics are known to be highly nonlinear, particularly when including the engine dynamics. However, in all longitudinal control architectures, the control is performed in two loops: an inner loop which compensates the nonlinear vehicle dynamics and achieves the desired algebraic acceleration by acting on the throttle or on the brakes, and an outer loop which is responsible for guaranteeing good tracking of the desired vehicle inter-distance. This decomposition is also called vehicle dependent and vehicle independent control parts. The vehicle-dependent loop is beyond the scope of this paper. It is assumed that it has already been designed and we are only interested here in the vehicle-independent part.

Longitudinal control includes the cruise control mode functionality which ensures the tracking of the desired cruise velocity v c . It switches automatically to the ACC mode when a preceding vehicle appears in the range of the host vehicle sensor with a slower speed.

Given the positions of the target and the following vehicles x t (t) and x f (t) respectively, the headway distance is s = x t (t)x f (t) and the relative speed is v r = ẋt (t)ẋ f (t). The headway distance dynamics can be represented as a double integrator driven by the difference between the leader vehicle acceleration and the follower vehicle acceleration, i.e.,

s = ẍl (t) -ẍ f (t) (1)

B. Spacing reference models

Several spacing strategies have been tested both for ACC systems and in the scope of the fully automated longitudinal control for automated highways. They include constant spacing, constant time headway and nonlinear time or space dependent headway.

1) Constant time headway: During the assistance intervention, the control law aims at tracking a constant time headway distance. This inter-distance policy is proportional to the following vehicle speed, v f (t), which leads to more spacing at higher speed. This is also very convenient since it can be implemented using only on board sensors. Desired spacing is thus given by s d (t) = s 0 + hv f (t), where s 0 is the minimum spacing at stopping and h is the time headway.

The spacing error, δ (t) = s(t)s d (t), represents the difference between the actual spacing and the desired one. The final state-space model is given by:

δ (t) = v r (t) -hu dr (t) vr (t) = -u dr (t) + γ l (t) (2)
where u dr is the acceleration/brake input set by the driver. When the assistance is active, it is assumed that an x-by-wire system allows canceling the driver inputs (acceleration and braking) and thus the acceleration/braking input is that of the controller u a . In that case, the model of eq. ( 2) reduces to:

ẋa (t) = Āx a (t) + Bu u a (t) + Bγ γ l (t) where (3)
Ā = 0 1 0 0 , Bu = -h -1 , Bγ = 0 1 .
The state vector is x a (t) (δ (t), v r (t)) T . The input γ l (t) = ẍl (t) is the leading vehicle acceleration which acts in this case as an unknown disturbance input. The whole state vector is used for control synthesis. The two variables are measured using on board sensors. In order to obtain a controller that is directly implementable on an Electronic Control Unit (ECU), the continuous-time state-space equation ( 3) is converted into a discrete-time model using a zero-order hold assumption on u a (t) and an exact discretization method with sample time T s . The sampling times t = kT s , where k represents the successive time steps belonging to I , the set of integers,

x(k + 1) = Ax(k)+B u u a (k) + B γ γ l (k) where (4) A = 1 T s 0 1 , B u = -hT s - T 2 s 2 -T s , B γ = T 2 s 2 T s . (5)
Vehicle longitudinal control needs to incorporate constraints on the absolute value of the jerk j f (t) of the host vehicle. In addition, an integral action in the open-loop ensures zero steady state error, thus the model of equation ( 4) is converted into an incremental input-output model (IIO). The new control input is in this case u(k) = δ u a (k) = u a (k)u a (k-1) and the IIO model is obtained by extending the state to x e (k) = (δ (k), v r (k), u a (k -1)) T . In addition, it is assumed that the leading vehicle acceleration is bounded, |γ l | ≤ γ max l , one defines a normalized disturbance input |w| ≤ 1 such that γ l = γ max l w and then the IIO model takes the form:

x e (k + 1) = A e x e (k) + B eu u(k) + B ew w(k) ( 6
)
where

A e =   1 T s (-hT s - T 2 s 2 ) 0 1 -T s 0 0 1   (7) B eu =   (-hT s - T 2 s 2 ) -T s 1   , B ew =   T 2 s 2 γ max l T s γ max l 0   (8)
In order to tackle with various traffic conditions, the headway time is assumed to be variable in the range [h mim , h max ]. Matrices in system [START_REF] Lingyun | A comprehensive review of the development of adaptive cruise control systems[END_REF] are thus functions of h, namely A e (h) and B eu (h). It could thus be converted into an LPV system with two vertices system obtained for the extremal values of h. One can write:

x e (k + 1) = 2 ∑ i=1 µ i [A ei x e (k)+B eui u(k)] +B ew w(k) (9)
where

[A e1 , B eu1 ] = [A e (h min ), B eu (h min )], [A e2 , B eu2 ] = [A e (h max ), B eu (h max )].
The coefficients µ 1 and µ 2 are positive and are verify equations (µ 1 + µ 2 = 1) and (µ 1 h min + µ 2 h max = h), which determine uniquely the parameters.

Our aim is now to synthesize an LPV state feedback of the form

u(k) = 2 ∑ j=1 µ j F j x e (k) (10)
By combining ( 9) and ( 10), the augmented closed-loop system is given by:

x e (k + 1) = 2 ∑ i, j=1
µ i µ j (A ei + B eui F j ) x e (k)+B ew w(k) (11) In the following, we denote Φ = ∑ 2 i, j=1 µ i µ j (A ei + B eui F j ), the closed loop system takes the form:

x e (k + 1) = Φx e (k) + B ew w(k) (12)

C. Invariant set under state feedback

Assume that there exists a quadratic function V (x e ) = x T e Px e , where P is a symmetric, positive definite matrix that satisfies, for all x e , w satisfying (12), w T w ≤ 1, V (x e ) ≥ 1, the condition [START_REF] Blanchini | Set-theoretic methods in control[END_REF]:

V (x e (k + 1)) ≤ V (x e (k)) (13)
Consider the reachable set Λ defined by:

Λ {x e (k ′ )| x e , w satisfying (12), x e (0) = 0, w T w ≤ 1, k ′ ≥ 0} (14)
The set ε P defined by:

ε P = {x e (k) ∈ R 3 |x e (k) T Px e (k) ≤ 1}, (15)
is an invariant set for system [START_REF] Naus | Design and implementation of parameterized adaptive cruise control: An explicit model predictive control approach[END_REF] with w ∈ R, w T w ≤ 1. This means that every trajectory that starts inside ε P remains inside it for k ≥ 0. Similarly, we will define the set ε ηP where η ∈ (0, 1] by

ε ηP = {x e (k) ∈ R 3 |x e (k) T ηPx e (k) ≤ 1} (16)
The existence of such a function V (x e ) means that the set ε P is an outer approximation of the reachable set Λ.

ε P is also an outer approximation of the reachable set

Λ * {x e (k ′ )| x e , w satisfying equation (12), x e (0) ∈ ε P , w T w ≤ 1, k ′ ≥ 0} (17)
In this section, the control law and the invariant set ε P are synthesized. This is achieved using BMI (Bilinear Matrix Inequalities) optimization method such that the system without the disturbance is asymptotically stable and at the same time, the reachable set for an initial state values inside the invariant set is contained in this invariant set.

D. Invariant set -quadratic boundedness

According to the previous considerations, the closed loop linear system x e (k + 1) = Φx e (k) + B ew w(k) is strictly quadratically bounded with a common Lyapunov matrix P > 0 for all allowable w(k) such that w(k) T w(k) ≤ 1, for t > 0, if x e (k) T Px e (k) > 1 implies (Φx e (k) + B ew w(k)) T P (Φx e (k) + B ew w(k)) < x T e (k)Px e (k), for any w T w ≤ 1.

The corresponding condition is obtained using the S-procedure and invoking the Schur complement [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. We want to impose the condition w(k) T w(k) ≤ 1 and

x T e (k)Px e (k) ≥ 1 implies V (x e (k + 1)) < V (x e (k)). A slightly more conservative condition that is more amenable to LMI formulation is to impose w(k) T w(k) ≤ x T e (k)Px e (k) implies V (x e (k + 1)) < V (x e (k)). This results in [4]: 2 ∑ i, j=1 µ i µ j ϒ i j > 0 (18)
where

ϒ i j = -Φ T i j PΦ i j + P -αP -Φ T i j PB ew -B T ew PΦ i j αI -B T ew PB ew (19)
and α > 0.

Using the Schur complement, and writing F j = Y j G -1 , where G = P -1 , equation ( 18) is equivalent to the inequality ∑ 2 i, j=1 µ i µ j Ῡi j > 0. Matrices Ῡi j are given by

Ῡi j =     G 0 αG GA T ei +Y T j B T eui 0 αI 0 B T ew αG 0 αG 0 A ei G + B eui Y j B ew 0 G     (20)
Finally, it is known that ∑ 2 i, j=1 µ i µ j Ῡi j ≥ 0 holds if there exist matrices Θ i j , i, j = 1, 2, Θ 11 = Θ T 11 and

Θ 12 = Θ T 21 such that            Ῡ11 > Θ 11 Ῡ22 > Θ 22 Ῡ12 + Ῡ21 ≥ Θ 12 + Θ 21 Θ 11 Θ 12 Θ 21 Θ 22 ≥ 0 (21)

E. Constraints on input and state components

In addition, it is possible to handle constraints on the control signal and the state:

-ū ≤ u(k) ≤ ū, -Ψ ≤ Ψx e (k + 1) ≤ Ψ, ∀t ≥ 0 (22)
where ū > 0 and Ψ := [ Ψ1 , Ψ2 , ..., Ψq ] T with Ψm > 0, m = 1, . . . , q, Ψ ∈ R q×3 and q is the number of imposed constraints. Notice that the bounds are provided separately on each state variables or a combination of them.

For a pre-specified scalar η ∈ (0, 1], the quadratic boundedness property ensures that if x e (0) ∈ ε ηP , then x e (k) ∈ ε ηP , ∀k ≥ 0, thus ∀w(k) such that w(k

) T w(k) ≤ 1 max k≥0 |u(k)| 2 = max k≥0 ∑ 2 j=1 µ j F j x e (k) 2 ≤ max k≥0 ∑ 2 j=1 µ j F j [ηP] -1/2 2 × [ηP] 1/2 [x e (k)] 2 ≤ ∑ 2 j=1 µ j F j [ηP] -1/2 2 ≤ ū2 (
ηG Y T j Y j ū2 ≥ 0 (24)
A similar procedure can be applied for the constraints on the state variables. Define ξ k as the k-th row of the q-ordered identity matrix. Then

max k≥0 |ξ k Ψx e (k + 1)| 2 = max k≥0 ξ k Ψ Φ B ew x e (k) w(k) 2 ≤ 2 ξ k Ψ Φ B ew ηP 0 0 I -1/2 2 (25)
If there exists a symmetric matrix Ξ such that

Ξ -2Ψ Φ B ew ηP 0 0 I -1 Φ T B T ew Ψ T ≥ 0 (26)
Ξ mm ≤ Ψ2 m , m = 1, ..., q, then Ψx e (k + 1) < Ψ, ∀k ≥ 0. By applying the Schur complement, and by pre-and postmultiplying the left-hand side of the obtained inequality with diag{G, I, I}, it is shown that (26) is equivalent to

       ∑ 2 i, j=1 µ i µ j   ηG * * 0 I * √ 2Ψ A ei G + B eui Y j √ 2ΨB ew Ξ   ≥ 0 Ξ mm ≤ Ψ2 m , m = 1, ..., q, ( 27
)

F. Inclusion of the normal zone

For the beginning, the qualitative description of the safety zones has been fitted in a mathematical form. The green zone corresponds to a predefined strip of length 2δ N centered on the reference distance s d (see Figs. 1(a) and 1(b)). In addition, the relative speed has been supposed to be bounded. These two bounded regions have been considered symmetrical with respect to the origin and denoted by |v r | ≤ v N r and |δ | ≤ δ N . The bound on the vehicle acceleration/brake is handled through the |u a | ≤ u N a . Putting together the above constraints, by definition, during "normal following" the state x e stays in a bounded space region:

|x ei | ≤ x N ei for i = 1, 2,
3, where x ei denotes the i-th component of the state vector x e which then belongs to the set

L(F N ) {x e ∈R 3 : | f N i x e | ≤ 1, i = 1, 2, 3}, (28)
where

F N ∈ R 3×3 , f N i represent the rows of F N , f N i,i = (x N ei ) -1 and f N i, j = 0 for i = j, i, j = 1, 2, 3. L(F N
) represents a polyhedron characterized by the diagonal matrix F N . The eight vertices of the normal zone

z N i = ±δ N , ±v N r , ±u N a T
have to be inside ε P , which leads to the LMI constraints:

1 (z N i ) T z N i G 0, i = 1, ...8 (29)

G. Avoiding entering the collision zone

Similarly, the condition that the vehicle does not enter the red zone (see Fig. 1(a)) can be written as |δ | ≤ δ M . Moreover, if the relative speed is limited, |v r | ≤ v M r , and the required longitudinal acceleration/deceleration is remaining within acceptable comfort values |u a | ≤ u M a , the three above constraints also define a parallelepiped denoted by L(F M ). A "safe following" is then defined by x e ∈ L(F M ), set which is greater than L(F N ). The region of the state space between these parallelepipeds corresponds to the orange zone, while the near region outside L(F M ) is the red zone. The matrix F M is given by:

F M = diag{ δ M -1 , v M r -1 , u M a -1 } = f M 1 f M 2 f M 3 , and L(F M ) {x e ∈ R 3 : | f M i x e | ≤ 1, i = 1, 2, 3}. (30)
The conditions ensuring that the invariant set ε P in inside this hypercube are [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]:

f M i G f M i T -1 < 0, i = 1, 2, 3. (31)

H. Optimization problem for control synthesis

We are now able to provide the BMI optimization problem with decision variables G, Y j , Θ i j , (i, j = 1, 2), Ξ and α, which ends with the control gain vectors F 1 and F 2 : 21), eq.( 24), LMIs f rom eq.( 27), eq.( 29), eq.( 31)

min trace(G) G ≻ 0, α > 0, eq.(
(32)
This BMI optimization problem is solved using the software developed by the authors of [START_REF] Kocvara | PENBMI Users Guide (Version 2.1), software manual[END_REF]. The trace of matrix G is minimized in order to reduce the excursions from the normal zone. It is also possible to maximize the trace of this matrix. This will make the assistance able to safely activate far from the normal zone. Notice that it is also possible to weight the magnitude of the maximum admissible disturbance which is also a design parameter.

IV. SIMULATION RESULTS

The assistance is now tested in a simulation environment which includes a vehicle dynamic model coupled to a driver model for car following, derived from the optimal velocity model of [START_REF] Bando | Analysis of Optimal Velocity Model with Explicit Delay[END_REF].

A. Numerical values

The "normal following" zone is defined with the following numerical values corresponding to an error of 10% of the nominal values for the speed v t f and inter-distance s t d = s 0 + hv t f , according to traffic situation: v N r = 0.1v t f , δ N = 0.1s t d . The "safe following" zone numerical values are chosen using the same idea but with different proportions: v M r = 0.5v t f , δ M = 0.7s t d . The optimization problem is solved while the BMI and LMI constraints are found feasible. The achieved numerical values are: α = 0.05, F 1 = 0.1839 0.1166 -0.0825 F 2 = 0.4398 -0.2786 -0.0995

G =   2.

 

The assistance is mainly designed in order to operate inside soft braking and acceleration modes. The brake and the acceleration bounds are chosen equal in absolute values to u max = 6m/s 2 , while the leading vehicle maximum brake magnitude is γ max l = 10m/s 2 . The solutions for the BMI optimization problems from eq. ( 32) lead to the invariant set shown in Fig. 2. This set includes the normal driving zone and is inside the limit zone.

B. Simulation results for braking and acceleration profiles

The distance at standing is set to s 0 = 5m and the time headway to h = 2 sec. Initial conditions are s(0) = 75m, v f (0) = 30m/s and v l (0) = v t f = 20m/s. Figure 3 illustrates the behavior of the assisted vehicle controlled by the control law. The test profile includes cruise control, hard stop and stop-and-go scenarios. The dash-dot lines in Fig. 3 correspond to the curves produced by the simulated motion of the lead vehicle. The solid line corresponds to the curves of the vehicle under the control assistance while the dash-dot line corresponds to the responses of the leader. The dashed curve in the first figure of Fig. 3 (a) corresponds to the reference inter-distance. When the follower vehicle comes near the leader vehicle, the velocity is adapted with comfortable deceleration and the vehicle is positioned to a safe distance. Then, at t = 25 sec the leader vehicle aims stopping with high braking value of 8.5m/s 2 , while the following vehicle comes to low speed of 5m/s and a corresponding headway distance of 15m. The critical distance s 0 = 5m is never reached. The braking is lower than 6m/s 2 . Thereafter, the leader vehicle is accelerated and decelerated (stop-andgo) with usual acceleration values but high jerk; however, the following vehicle is maintained to a safe distance, and with bounded jerk (< 3m/s 3 ). One can note from Fig. 2 that trajectory in the phase-plane starts at the limits of the ellipsoid and thus remains all the time inside.

In practice, the estimation of the leader vehicle acceleration could be affected by noise or a bias. Thus, the previous maneuver is simulated with non perfect measurement of the leader vehicle assuming a zero mean Gaussian additive noise signal with variance equal to 0.1 and a bias equal to 0.1m/s 2 . Results for the controlled vehicle are visible in Fig. 3. All the responses are smooth, the control law achieves a high quality filtering of the noise.

C. Testing for different headway time values

Suppose now that starting with a time headway of 1.5 sec, the infrastructure manager decides to linearly increase it until 2.5 sec between [20, 35] sec and hold it for a time duration of 20 sec, Thus it is linearly decreased again towards 1.5 sec between [55, 70] sec and then maintained to this value. The reference distance which was initially of 35m is now needed to increase to 55m. Figure 1 shows the reference distance profile and the corresponding experienced inter-distance. One can notice that the assistance adapts for this profile as the two plots are very close. The adaptation is smooth since it is performed with a speed variation of less than 1.2m/s.

V. CONCLUSION

This paper presents the design and the simulation test of a vehicle following assistance for variable headway time handling and. In fact, a control law based on linking Lyapunov functions and invariant sets for multi-model systems with LMI and BMI optimization is proposed. This approach ensures four important features: asymptotic convergence of the vehicle trajectory to the reference headway distance with variable headway time, guarantee that the vehicle remains inside a "safe following" zone, and a bounded control input for all bounded leading vehicle acceleration. Simulation tests show that entering again in the collision red zone is avoided once the invariance set is reached. It has been verified that safety and comfort for the entire range of speed are achievable. The field test of the proposed strategy on prototype vehicles including implementation in platooning strategy is the object of future work.

Fig. 1 .

 1 Fig. 1. (a) Mathematical representation of vehicle following zones, (b) Vehicle following zones.

  23)By applying the Schur complement, and pre-and postmultiplying the left-hand side of the inequality the obtained matrix with diag{G, I}, the condition (23) is satisfied if

Fig. 2 .

 2 Fig. 2. Following vehicle trajectories in the (Relative speed / Inter-distance error), phase plane.

Fig. 3 .

 3 Fig. 3. Follower response to leader profile in nominal conditions.

Fig. 4 .

 4 Fig. 4. Constant speed: Inter-distance, velocities, acceleration and jerk for variable headway time.

  696 × 10 2 1.687 × 10 2 6.673 × 10 2 1.687 × 10 2 1.297 × 10 2 3.975 × 10 2 6.673 × 10 2 3.975 × 10 2 4.516 × 10 3


Mengchu Zhou
email: zhou@njit.edu

P Wu

F Chu

John C Hartmann

Bi-Objective Scheduling of Fire Engines for Fighting Forest Fires: New Optimization Approaches

Peng Wu, Feng Chu, Senior Member, IEEE,A d aC h e ,Senior Member, IEEE, and MengChu Zhou, Fellow, IEEE Abstract-It is challenging to perform emergency scheduling for fighting forest fires subject to limited rescue resources (i.e., vehicles with fire engines), since extinguishing each fire point should take into account multiple factors, such as the actual fire spreading speed, distance from fire engine depot to fire points, fire-fighting speed of fire engines, and the number of dispatched vehicles. This paper investigates a bi-objective rescue vehicle scheduling problem for multi-point forest fires, which aims to optimally dispatch a limited number of fire engines to extinguish fires. The objectives are to minimize the total fireextinguishing time and the number of dispatched fire engines. For this problem, we first develop an integer program that is an improved and simplified version of an existing one. After exploring some properties of the problem, we develop an exact dynamic programming algorithm and a fast greedy heuristic method. Computational results for a real-life instance, and benchmark and large-size randomly generated instances confirm the effectiveness and efficiency of the proposed model and algorithms. Besides, a bi-objective integer program is developed to address the multi-depot fire engine scheduling issue. Index Terms-Forest fires, emergency scheduling, optimization, multi-objective, Pareto front.

I. INTRODUCTION

I N RECENT years, natural disasters and man-made catastrophic events happen frequently, which have caused tremendous harm to human beings. Effective and efficient emergency management is crucial in ensuring timely rescue and relief operations. Currently, the main studies about emergency management focus on emergency response management [START_REF] Ardekani | Logistics problems in the aftermath of the 1985 mexico city earthquake[END_REF]- [START_REF] Wang | An integrated road construction and resource planning approach to the evacuation of victims from single source to multiple destinations[END_REF], relief distribution [START_REF] Knott | Vehicle scheduling for emergency relief management: A knowledge-based approach[END_REF]- [START_REF] Misra | Priority-based time-slot allocation in wireless body area networks during medical emergency situations: An evolutionary game-theoretic perspective[END_REF], and emergency logistics management [START_REF] Ardekani | Transportation operations following the 1989 loma prieta earthquake[END_REF]- [START_REF] Saber | Resource scheduling under uncertainty in a smart grid with renewables and plug-in vehicles[END_REF]. To deal with the related problems, a number of modelling approaches, e.g., integer programming [START_REF] Sheu | An emergency logistics distribution approach for quick response to urgent relief demand in disasters[END_REF], [START_REF] Hu | A container multimodal transportation scheduling approach based on immune affinity model for emergency relief[END_REF], fuzzy optimization [START_REF] Zheng | Emergency transportation planning in disaster relief supply chain management: A cooperative fuzzy optimization approach[END_REF], and dynamic programming [START_REF] Yi | A dynamic logistics coordination model for evacuation and support in disaster response activities[END_REF], and solution techniques, e.g., simulation [START_REF] Haghani | Simulation model for real-time emergency vehicle dispatching and routing[END_REF], [START_REF] Fu | Optimization of evacuation traffic management with intersection control constraints[END_REF] and meta-heuristic [START_REF] Yi | Ant colony optimization for disaster relief operations[END_REF], have been proposed. Recent surveys related to various emergency operations can refer to [START_REF] Simpson | Fifty years of operational research and emergency response[END_REF]- [START_REF] Caunhye | Optimization models in emergency logistics: A literature review[END_REF].

Many highlighted works have adequately addressed emergency planning and scheduling problems under various disasters and emergency events such as earthquakes, hurricane, tsunamis, and floods. However, limited focus has been on the forest fires. Notably, with the warming of global climate, forest fires around the world become increasingly frequent. As reported in [START_REF] Wei | The usage state and the prospect of forest fire extinguish equipments[END_REF], the forest fires occur over 22 million times each year in the world and burn more than 6.4 million hectares of forest area, accounting for more than 0.23% of the world's forest coverage. In China, the forest area engulfed by forest fire each year is about 1.1 million hectares, accounting for about 0.8% to 0.9% of the whole national forest area. Moreover, many forest fires involve multiple fire points simultaneously. For example, the forest fires occurred in Nanwenghe region located in Mt. Daxing'anling on March 19, 2003 and Huzhong region located in Mt. Daxing'anling on June 29, 2010 involve four and seven fire points, respectively. Forest fire has been recognized as one of the eight major natural disasters [START_REF] Wei | The usage state and the prospect of forest fire extinguish equipments[END_REF]. The basic forest fire-fighting workflow is illustrated in Fig. 1 [START_REF] Zhang | The Standardized Operational Workflow of Forest Fire Fighting in Daxing'Anling Area[END_REF]. Once forest fires happen, the fire information is collected and the fire level is determined. Then, the emergency plan is launched. Next, the fire-fighting plan is executed, which usually involves five main groups (i.e., organization and command, emergency communication, fire engines scheduling, logistics support and medical service, and propaganda), followed by fighting fires and reporting the results. Finally, the firegrounds are cleaned and guarded and the post-disaster disposal is conducted.

In this study, we focus on the fire engine scheduling that is finished by the fire engines scheduling group as indicated by dashed line rectangle in Fig. 1. For each fire engine dispatched, it is usually assigned with an identical quantity of firemen and fire-fighting equipment [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. The most critical question facing the fire engine scheduling group lies in how to dispatch fire engines to extinguish the fire points in the most efficient manner. Therefore, it is necessary to develop effective scheduling approaches to achieve the highest firefighting efficiency while reducing or limiting the expenditures.

Different from the other natural disasters, once forest fires happen, different fire points may spread at various speeds that are usually determined by natural conditions, e.g., wind force, terrain slope, and fuel types [START_REF] Wang | The measurement method of the wildfire initial spread rate[END_REF]. For example, after the forest fires occurred in Mt. Daxing'anling, China in 2010, its seven fire points spread at different speeds varying between 2.2 m/min and 6.98 m/min [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. Different forest fire starting time would affect the actual fire engine scheduling. For example, if the starting time intervals between different fire points are large, decision-makers can address them one by one. However, if the starting time intervals between different fire points are short, the fire engine scheduling should consider multiple fire points simultaneously. In addition, extinguishing each fire point should take into account multiple other factors such as the fire spreading speed, distance from fire engines depot to fire points, dispatched fire engine count and their fire-fighting speed. Consequently, due to these special characteristics, the existing methodologies for emergency issues like earthquakes cannot be directly applied to the emergency response scheduling to tackle forest fires. Furthermore, we note that most of the existing studies mainly focused on minimizing the emergency response time or response costs. However, to make a wellinformed decision, decision-makers may desire to know the service performance corresponding to different resource usage, i.e., the trade-off between the two aspects. To simultaneously optimize the emergency response service and resource usage can help decision-makers make a good trade-off between them. For example, they usually desire to minimize the total fire-extinguishing time and expenditures for dealing with the forest fires at the same time. Besides, a number of studies, e.g., [START_REF] Haghani | Formulation and solution of a multi-commodity, multimodal network flow model for disaster relief operations[END_REF], [START_REF] Haghani | Decision support system for snow emergency vehicle routing: Algorithms and application[END_REF], [START_REF] Haghani | Snow emergency vehicle routing with route continuity constraints[END_REF], have also pointed out that minimizing emergency response costs is as important as minimizing the emergency response time. Consequently, multi-objective optimization is required to support the desired decisionmaking for processing emergency issues.

To the best of our knowledge, no studies in the literature have addressed the emergency planning and scheduling for the forest fires except the recent work [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. Tian et al. pioneered in a bi-objective optimization model to simultaneously optimize the total fire-extinguishing time and the total number of fire engines dispatched [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. To solve the proposed model, they developed a bi-objective hybrid meta-heuristic algorithm to obtain a set of Pareto solutions. Instances with up to 50 fire points and 250 fire engines were solved. However, their study has three deficiencies. Firstly, their obtained solutions are not guaranteed to be the best (Pareto optimal) due to the nature of meta-heuristic algorithms. Secondly, the computational time of their method is relatively long to provide satisfactory solutions, especially for large-size problems, which may greatly increase the total time to extinguish fires. Thirdly, the work [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF] does not address the multi-depot fire engine scheduling issue. To overcome the above deficiencies, this work intends to develop effective and efficient approaches for solving the problem and address the multi-depot fire engine scheduling. The main contributions of this paper are as follows:

1) We propose an improved and simplified integer program with fewer variables and constraints. Comparison results show that the improved model significantly outperforms the existing one. Moreover, we extend the proposed model to the multi-depot case. 2) We prove that the optimal fire-extinguishing time is a strictly decreasing function of the number of dispatched fire engines. Then, we derive the lower and upper bounds on the number of dispatched fire engines and show that the number of non-dominated points over the Pareto front is their linear function. 3) Based on the derived problem properties, we develop an exact dynamic programming algorithm and a fast greedy heuristic method for solving the problem. Computational results for a real-life case, benchmark and randomly generated large-size instances confirm the effectiveness and efficiency of the proposed algorithms. Besides, the ε-constraint method is adapted to obtain the Pareto front of the multi-depot fire engine scheduling problem. The remainder of the paper proceeds as follows. In Section II, the related literature is reviewed. The problem description is recalled and an improved bi-objective integer program is provided in Section III. Section IV analyzes the properties of the problem. Based on them, an exact dynamic programming algorithm and a greedy heuristic method are developed in Section V. Section VI addresses the fire engine scheduling problem with multiple depots. In Section VII, computational experiments are conducted. Conclusions and future research issues are discussed in Section VIII.

II. LITERATURE REVIEW

Emergency rescue and relief operations are mainly dependent on surface transportation. As a consequence, well-planned roadway networks and efficient rescue vehicle scheduling after a natural disaster or catastrophic event takes place is critical to ensure response efficiency and reduce costs. Over the past decades, a number of studies concentrated on emergency planning and scheduling for various natural disasters and catastrophic events. For instance, Haghani and Oh [START_REF] Haghani | Formulation and solution of a multi-commodity, multimodal network flow model for disaster relief operations[END_REF] presented a multi-commodity, multi-modal network model for determining the rescue vehicle routing and scheduling when facing natural disasters (e.g., earthquakes and hurricanes) to minimize the cost. Two heuristics were developed to solve it. Later, Özdamar et al. [START_REF] Özdamar | Emergency logistics planning in natural disasters[END_REF] presented a rescue vehicle routing model for optimally dispatching rescue commodities (e.g., medical materials and personnel) in natural disasters, e.g., earthquakes, to minimize the unsatisfied demands. A Lagrangean relaxation heuristic was proposed to yield near-optimal solutions. Campbell et al. [START_REF] Campbell | Routing for relief efforts[END_REF] formulated a rescue vehicle routing problem in aftermath of a natural disaster (e.g., earthquake, hurricane, and tsunami) as variants of the traveling salesman problem (TSP) and vehicle routing problem (VRP). Effective heuristics were proposed as well. Yuan and Wang [START_REF] Yuan | Path selection model and algorithm for emergency logistics management[END_REF] constructed a path selection model to minimize the total rescue time in emergency logistics management and proposed a modified Dijkstra algorithm for it. Zhang et al. [START_REF] Zhang | Route selection for emergency logistics management: A bio-inspired algorithm[END_REF] developed a bio-inspired algorithm to address the path selection in an emergency logistics network. Moreover, Wex et al. [START_REF] Wex | Emergency response in natural disaster management: Allocation and scheduling of rescue units[END_REF] constructed a rescue unit allocation model to process incidents caused by natural disasters, e.g., earthquakes, tsunamis and floods with the objective of minimizing the total completion time. Several heuristics were suggested to solve it.

Particularly, to effectively respond to earthquakes, Hu et al. [START_REF] Hu | A modified particle swarm optimization algorithm for optimal allocation of earthquake emergency shelters[END_REF] presented an optimal earthquake emergency shelter allocation model to optimize the total operating cost while meeting all victims' needs. A particle swarm optimization based metaheuristic was suggested to tackle the model. Yan et al. [START_REF] Yan | Optimal scheduling of logistical support for an emergency roadway repair work schedule[END_REF] developed a logistical support scheduling model for damaged roadway repair work after an earthquake to minimize the total operating cost. Later, Yan et al. [START_REF] Yan | Optimal scheduling for highway emergency repairs under large-scale supply-demand perturbations[END_REF] established an emergency scheduling model for highway emergency repairs under demand-side perturbations to minimize the overall completion time of repairing all highway segments. An ant colony system based hybrid meta-heuristic algorithm was developed to derive near-optimal solutions.

Besides, Haghani et al. [START_REF] Haghani | Simulation model for real-time emergency vehicle dispatching and routing[END_REF] developed a simulation model for optimal emergency vehicle dispatching and routing during accidents to minimize the average response time. Later, the work [START_REF] Yang | Online dispatching and routing model for emergency vehicles with area coverage constraints[END_REF] presented an online dispatching and routing model for optimally assigning emergency response vehicles to incidents as well as determining the routes to minimize the total emergency response time. Recently, Park et al. [START_REF] Park | A stochastic emergency response location model considering secondary incidents on freeways[END_REF] constructed a stochastic emergency response location-allocation model for processing traffic incidents on highways with the objective of minimizing the expected delay of all scenarios. A heuristic was developed to solve it. Haghani and Qiao [START_REF] Haghani | Decision support system for snow emergency vehicle routing: Algorithms and application[END_REF], [START_REF] Haghani | Snow emergency vehicle routing with route continuity constraints[END_REF] developed salting truck routing models for dealing with snow emergencies to minimize the number of trucks used and deadhead time. Constructive heuristics were proposed to obtain satisfactory solutions.

From the existing research results, it can be concluded that the prior studies on emergency planning and scheduling mainly focus on processing emergency issues including earthquakes, hurricane, tsunamis, floods, snow emergencies, and accidents. As previously mentioned, the forest fire has its own special characteristics. The recent work [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF] pioneered in the first multi-objective model and hybrid meta-heuristic for addressing the fire engine scheduling problem. However, the exisiting model is very hard to solve and the existing algorithm is time-consuming and may fail to find the Pareto front. In this research, we improve the existing research result [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF] in a number of ways by improving the formulation, developing new efficient exact and heuristic algorithms, and extending the studied one-depot problem to a multi-depot case.

III. PROBLEM DESCRIPTION AND FORMULATION

The emergency scheduling problem for forest fires can be defined over an undirected graph G = (V ′ , E),w h e r eV ′ = {0}∪V ={1,...,n} (resp. E) is the set of vertices (resp. the set of edges). Without loss of generality, vertex 0 represents the starting point that hosts all fire engines and vertex i ∈ V denotes the i -th forest fire point to be extinguished.

As is widely acknowledged, once a major fire disaster with multiple fire points happens simultaneously in a region, to extinguish all the fire points as quickly as possible is desired for each decision-maker. The whole rescue time for fire point i ∈ V is composed of the travel time of fire engines, denoted by T Ai , and extinguishing time, denoted by t Ei . T Ai is assumed to be calculated by the distance between the starting point 0 and fire point i , denoted by d 0i , and the average driving speed of fire engines, denoted by v 0i . T Ai = d 0i /v 0i in this paper. However, it is much more complex and difficult to calculate t Ei , since it depends on multiple factors, such as the fire spread speed v si , travel time of fire engines T Ai , fire-fighting speed of fire engines v f and dispatched vehicle count y i . In the paper, we compute t Ei using the following formula as proposed by [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]:

t Ei = v si T Ai /(y i v f -v si ), ∀i ∈ V
It can be seen that the fire lasting time is influenced by vehicle scheduling. Conversely, it would influence vehicle scheduling, since one objective is to minimize the total fire lasting time.

Remark 1: For any fire point i ∈ V ,thev alueofy i v f -v si should be greater than 0; otherwise this fire point cannot be extinguished. In other words, y i >v si /v f should be ensured in processing fires to derive a feasible scheme.

The fire spread speed for any given fire point is determined by the following fire spread model taking into account natural factors, e.g., wind force, terrain slope, and fuel types [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF].

v s = v 0 λ s λ ϕ λ w = v 0 λ s λ ϕ e 0.1783v w
where v s is the fire spread speed; v 0 = αT + βw + γ denotes the initial spread speed; λ s ,λ ϕ and λ w are the correction factors of fuel types, terrain slope and wind force, respectively; T and w represent the temperature and wind force, respectively; α, β and γ are the factors related with the terrain, respectively. According to the above formula for computing t Ei , the fire spread speeds of different fire points, denoted by v si , can greatly differ among each other due to different natural conditions. For more details on the above formulas, see [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. Since the total number of fire engines is limited, it can be imagined that if one fire point is assigned with more fire engines, the available fire engines to be dispatched to the rest becomes fewer. However, different fire points are independent such that the spread of different fire points would not influence each other. As desired by decision-makers, to minimize the total fire-extinguishing time is the primary goal. Meanwhile, the dispatched fire engine count is desired to be optimized as the second objective. This study aims to optimize the above two objectives simultaneously via an optimal emergency schedule such that an appropriate number of fire engines are dispatched to different fire points.

A. Problem Formulation

With the notations presented in Table I, we first recall the following integer program presented in [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF].

P 0 : f 1 : min n i=1 t Ei
(1)

f 2 : min n i=1 m j =1 x ij (2) s.t. t Ei = v si T Ai /( m j =1 x ij v f -v si ), ∀i ∈ V (3) LB ≤ n i=1 m j =1 x ij ≤ m (4) L i ≤ m j =1 x ij ≤ U i , ∀i ∈ V (5)
x ij ∈{0, 1}, j = 1, 2,...,m, ∀i ∈ V (6)

t Ei ≥ 0, ∀i ∈ V (7)
Objective ( 1) is to minimize the total fire-extinguishing time.

Objective ( 2) is to minimize the number of dispatched rescue vehicles. Constraint (3) defines the fire-extinguishing time of fire point i ∈ V . Constraint (4) aims to guarantee that the total number of fire engines sent to all fire points cannot be greater than the number of available fire engines m and less than LB vehicles such that all fire points can be extinguished. Constraint [START_REF] Wang | An integrated road construction and resource planning approach to the evacuation of victims from single source to multiple destinations[END_REF] restricts the number of fire engines dispatched to each fire point. Constraints ( 6) and ( 7) impose restrictions on decision variables. It can be observed that model P 0 employs mn integer variables and 3n + 2 constraints.

B. An Improved Formulation

In this section, we improve the above formulation. Firstly, we can directly employ integer variable y i to replace binary variable x ij that can reduce (m -1)n variables, since the fire engines are supposed to be indentical [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. Secondly, we observe that parameter LB in P 0 is redundant because it can be implicitly ensured by the left side of ( 5). Thus we can remove such bound constraint of (4). Thirdly, the right side of (5) can also be removed since it can be ensured through limiting the total number of fire engines as done in (4). With these above observations, the following equivalent improved and simplified formulation can be derived.

P : f 1 : min n i=1 t Ei f 2 : min n i=1 y i (8) s.t. t Ei = v si T Ai /(y i v f -v si ), ∀i ∈ V (9) n i=1 y i ≤ m ( 10
)
L i ≤ y i , ∀i ∈ V (11) t Ei ≥ 0, y i ∈ IN + , ∀i ∈ V ( 12
)
It can be seen that model P employs n integer variables and 2n +1 constraints. Compared with model P 0 , model P reduces (m -1)n integer variables and n + 1 constraints. Comparison results reported later show that P is much more efficient to solve than P 0 .

We note that both P 0 and P involve integer variables and are non-linear due to the non-linear constraints (3) and ( 9). It is difficult to solve them by commercial optimization software such as LINGO to be reported later. In order to solve P, its property analysis is performed in the following section.

IV. PROPERTIES OF THE MODEL

As a schedule with fewer rescue vehicles dispatched to extinguish fires may lead to a longer fire-extinguishing time, the two objectives of the proposed model (i.e., fireextinguishing time and number of dispatched fire engines) are conflicting. As such, there does not exist an unique solution to simultaneously optimize the above two objectives but a set of Pareto-optimal solutions. For simplicity, let Y denote the feasible solution set of P. Before proceeding, we first give the following definitions related to multi-objective optimization [START_REF] Wu | An improved exact ε-constraint and cut-and-solve combined method for biobjective robust lane reservation[END_REF].

Definition 1: A schedule y * ∈ Y is Pareto-optimal if and only if there exist no schedule y ∈ Y, such that f 1 (y) ≤ f 1 (y * ) and f 2 (y) ≤ f 2 (y * ) where at least one inequality is strict.

Definition 2: A solution (point) ( f 1 (y), f 2 (y)), y ∈ Y,i n the objective space is non-dominated if y is Pareto-optimal.

The image of all non-dominated points is called the nondominated set or the Pareto f ront. To exactly solve a multiobjective optimization problem means finding its Pareto front.

To derive the relation between the two objectives of model P, we first define the feasible (expected) number of dispatched fire engines as a parameter R satisfying n i=1 L i ≤ R ≤ m where m and L i is the total available fire engine count and the minimum number of fire engines needed to fight fire point i , respectively. Then, model P can be transformed into P R for a given R.

P R : min i∈V t Ei s.t. i∈V y i ≤ R ( 13
)
and constraints ( 9), ( 11) and ( 12) P R is to minimize the total fire-extinguishing time for a given number of rescue vehicles R. For the convenience of analysis, let

y * R = (y * 1,R , y * 2,R ,...,y * n,R ), t * Ei ,a n dT * R = n i=1 t *
Ei denote an optimal number of fire engines dispatched to fire points i = 1, 2,...,n, the time of extinguishing fire point i , and the optimal objective value, respectively. For any problem P R , we have the following results.

Proposition 1: In an optimal solution of P R , the number of dispatched fire engines is exactly equal to R, i.e., 9), for any i ∈ V ,if△ i > 0, then the value of t Ei in the new schedule must be less than that of the original one. This means that we have found a new schedule with less fire-extinguishing time and R rescue vehicles are dispatched to fight fires. This contradicts the fact that y * R is an optimal schedule. Consequently, n i=1 y * i,R = R always holds. Proposition 2: In an optimal solution, the fire-extinguishing time T * R is a strictly decreasing function of the number of dispatched fire engines R, i.e., for any two values

n
R 1 and R 2 , if R 1 < R 2 ,t h e nT * R 1 > T * R 2 always holds. Proof: We firstly show that if R 1 < R 2 ,t h e nT * R 1 ≥ T * R 2
always holds. For any given two values of R, i.e., R 1 and R 2 ,i fR 1 < R 2 , then constraint (13) of P R 1 is tighter than that of P R 2 and all the other constraints of both problems are the same. Thus,

T * R 1 is not less than T * R 2 , i.e., T * R 1 ≥ T * R 2 . Then, we show T * R 1 = T * R 2 is impossible by contradiction. Assume that T * R 1 = T * R 2 .
Then, the optimal schedule of P R 1 is also an optimal one of P R 2 . By Proposition 1, we have

n i=1 y * i,R 1 = R 1 = R 2 , which contradicts the fact that R 1 < R 2 . This indicates that it is impossible T * R 1 = T * R 2 . Consequently, if R 1 < R 2 ,t h e nT * R 1 > T * R 2 always holds.
By Propositions 1 and 2, we have the following corollary. Corollary 1: A Pareto-optimal solution can always be obtained by exactly solving P R for a given R.

Furthermore, as n i=1 L i ≤ R ≤ m, the lower and upper bounds on the number of dispatched fire engines (i.e., f 2 ) over the Pareto front are n i=1 L i and m, respectively. We can also have the following theorem.

Theorem 1: There exist m -n i=1 L i + 1 non-dominated points in the Pareto front of model P and they can be obtained by separately solving P R with R being each value in the set ={ n i=1 L i , n i=1 L i + 1,...,m}. Proof: By Propositions 1 and 2 and the fact that the value of f 2 must be an integer.

In the following, new optimization approaches based on the derived properties are developed for the considered problem.

V. S OLUTION APPROACH

As indicated by Theorem 1, the key to solving P becomes how to efficiently solve P R .S i n c eP R is a non-linear integer program because of constraint ( 9), it is difficult to be tackled by a commercial solver such as LINGO. Therefore, it is necessary to develop solution approaches for it. In what follows, we design an exact dynamic programming algorithm and a fast greedy heuristic method for solving P R .F o rt h e former, the Pareto front is ensured at the expense of computational time, while for the latter, an approximated Pareto front is derived but the computational time can be greatly saved.

A. Dynamic Programming Algorithm

To solve problem P R exactly, we develop in this subsection a dynamic programming (DP) algorithm based on the properties of the problem.

From P R , we can find that for each fire point i , a minimum number of rescue vehicles L i must be dispatched such that it can be extinguished (see constraint [START_REF] Misra | Priority-based time-slot allocation in wireless body area networks during medical emergency situations: An evolutionary game-theoretic perspective[END_REF]). Consequently, for P R , the key problem becomes how to dispatch the remaining R -n i L i vehicles such that the total fire-extinguishing time is minimized. Before proceeding, we first define a new variable as follows.

z i : additional number of fire engines except a minimum value L i dispatched to fire point i , i ∈ V . Note that y i = z i + L i . With the new definition and Proposition 1, P R can be rewritten as follows:

min n i=1 v si T Ai /((z i + L i )v f -v si ) s.t. n i=1 (z i + L i ) = R ( 14
)
z i ≥ 0, ∀i ∈ V ( 15
)
Let

T i = v si T Ai /(L i v f -v si )
, which denotes the fireextinguishing time of fire point i with L i fire engines. Then, the above program is equivalent to the following one:

P ′ R : max n i=1 (T i -v si T Ai /((z i + L i )v f -v si )) (16) s.t. n i=1 z i = R - n i=1 L i z i ≥ 0, ∀i ∈ V ( 17
)
To solve P ′ R , a DP algorithm is proposed as follows. For the convenience of analysis, let r i (z i ) = T i -v si T Ai / ((z i + L i )v f -v si ); and for fire point i ,l e tF i (z) denote the maximum objective value among all vehicle schedules that contain fire point set {i, i + 1,...,n},g i v e nz fire engines to be dispatched. By using the above notations, the optimal solution of P ′ R can be computed in a recursive way

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ F n (z) = r n (z) F i (z) = max 0≤z i ≤z {r i (z i ) + F i+1 (z -z i )}, ∀ 1 ≤ i < n (18)
The optimal objective value of P R can be obtained by

F 1 (R -n i=1 L i ) in O(nm 2
) time in the worst case. Thus, the optimal objective value of P R is n i=1 T i -F 1 (Rn i=1 L i ). By Theorem 1, there are m-n i=1 (at most m) single-objective problems P R 's need to be solved, thus the whole solution of model P with the DP algorithm requires O(nm 3 ) time in the worst case, which is pseudo polynomial.

B. Greedy Heuristic Algorithm

In order to more efficiently solve P, especially for large-size problems, we also propose a greedy heuristic method for P R whose core idea is to first ensure the minimum requirements of fire engines for each fire point and then to assign remaining available vehicles one by one to the fire point with the largest reduction in the total fire-extinguishing time. Its procedure is outlined in Algorithm 1.

Algorithm 1 A Greedy Heuristic Algorithm for

P R 1: Let y i = L i , ∀i ∈ V and obtain y = (L 1 , L 2 ,...,L n ) 2:
Calculate f 1 (y) and let R v be R -n i=1 L i , which denotes the remaining fire engine count expect the minimum requirement.

3: while (R v > 0) 4: Let Rt i = 0, ∀i ∈ V 5: for i = 1; i ≤ n do 6:
Let y i = y i +1 and obtain solution y ′ = (y 1 , y 2 ,...,y n )

7: Calculate f 1 (y ′ ) and let Rt i = f 1 (y) -f 1 (y ′ ) 8:
y i = y i -1 9: end for 10: Select the fire point i =a r gm a xRt i 11: y i = y i + 1a n dR v = R v -1 12: end while 13: Output the solution y = (y 1 , y 2 ,...,y n ) and the corresponding objective value f 1 (y).

In Algorithm 1, step 1 is to assign the minimum number of fire engines to each fire point; step 2 calculates the corresponding fire-extinguishing time and the rest fire engine count; steps 3-12 are to assign the remaining fire engines one by one where steps 5-9 are to determine which fire point to be assigned at each iteration; and step 13 outputs the solution and its corresponding objective value. Clearly, the most time-consuming parts of Algorithm 1 are steps 3-12. Since 0 ≤ R l ≤ mn, then at most mn additional rescue vehicles are to be assigned. Then, for each assignment (i.e., steps 5-9), we need to check at most n fire points. Thus, Algorithm 1 requires O(nm) time in the worst case.

As Algorithm 1 cannot guarantee the solution optimality, an approximated Pareto front of P is found by solving all P R 's using Algorithm 1. Similarly, since at most m problems need to be solved, the whole solution of model P using Algorithm 1 for each P R requires O(nm 2 ) time in the worst case, which enjoys lower time complexity than the DP algorithm.

VI. EXTENSION TO MULTI-DEPOT CASE

In practical situations, it might be possible that a single fire-fighting service point cannot meet all the fire-fighting needs, i.e., the capability of the fire-fighting service station is not enough. This motivates us to deploy multiple service points to collaboratively extinguish forest fires. This subsection expands model P to a multi-depot case, thus formulating a multi-depot fire engine scheduling problem for the first time.

To formulate this problem, one of the key problems is to calculate the fire-extinguishing time. For the sake of convenience, let K , m k , T Aik and y ik denote the set of all service points, the number of fire engines involved in the k-th depot, the travel time of fire engines from the k-th service point (depot) to fire point i , and the number of fire engines dispatched to fire point i from the k-th service point, respectively. Note that y ik is the decision variable. Similarly, T Aik can be calculated as d ki /v ki ,w h e r ed ki and v ki are the distance between depot k and fire point i and the average driving speed of fire engines, respectively [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. Let t Ri denote the whole rescue time for fire point i ∈ V . Then, according to [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF], the mathematic relation among v si , v f and T Ri can be described as follows:

t Ri • v si = k∈K (t Ri -T Aik ) • v f y ik , ∀i ∈ V ,
where the left (resp. right) side represents the fire-fighting (resp. fire spread) area of the i -th fire point and t Ri -T Aik denotes the fire-fighting time of the fire engines from the k-th depot. Thus, the whole rescue time for fire point i is:

t Ri = k∈K v f T Aik y ik k∈K v f y ik -v si , ∀i ∈ V.
Remark 2: For any fire point i ∈ V , k∈K y ik > v si /v f should be ensured to derive a feasible scheme.

Thus, the multi-depot fire engine scheduling problem can be formulated as the following program P M .

P M : f 1 : min i∈V t Ri (19) f 2 : min i∈V k∈K y ik (20) s.t. t Ri = k∈K T Aik v f y ik k∈K v f y ik -v si , ∀i ∈ V (21) n i=1 y ik ≤ m k , ∀k ∈ K (22) k∈K y ik ≥ L i , ∀i ∈ V ( 23
)
t Ri ≥ 0, y ik ∈{0}∪IN + , ∀i ∈ V , ∀k ∈ K ( 24
)
where constraint [START_REF] Fu | Optimization of evacuation traffic management with intersection control constraints[END_REF] ensures that the dispatched fire engines from the k-th depot should not exceed its capacity and constraint [START_REF] Yi | Ant colony optimization for disaster relief operations[END_REF] ensures that the number of fire engines dispatched to fire point i should exceed the given lower bound value L i ≥⌈ v si /v f ⌉ b a s e do nR e m a r k2 .A sP M has its own specialities, the proposed DP and heuristic for P cannot be used to solve it despite our best efforts. To solve P M ,t h eεconstraint method [START_REF] Haimes | On a bicriterion formulation of the problems of integrated system identification and system optimization[END_REF] is adapted to solve it. Its basic idea is to transform the bi-objective model to a single-objective one by optimizing a preferred objective while restricting the other to a given bound ε. For more details on the ε-constraint method, please refer to [START_REF] Wu | An improved exact ε-constraint and cut-and-solve combined method for biobjective robust lane reservation[END_REF] and [START_REF] Haimes | On a bicriterion formulation of the problems of integrated system identification and system optimization[END_REF]. For P M , we select the first objective f 1 as the preferred one. Then, P M can be transformed to P M (ε) shown as follows:

P M (ε) : min i∈V t Ri s.t. i∈V k∈K y ik ≤ ε and constraints (21) -(24) (25)
where ε is an upper bound of f 2 . It is not hard to find that the range of ε would be [ n i=1 L i , k∈K m k ].G i v e na value of ε, one non-dominated point may be found by exactly solving P M (ε). The Pareto front can be obtained by varying the value of ε within its range. Its procedure can be described in Algorithm 2.

Algorithm 2

The ε-Constraint Method for P M 1: Transform P M to P M (ε) 2: Calculate the range of f 2 , i.e., [ n i=1 L i , k∈K m k ]. 3: for ε ← n i=1 L i to k∈K m k do 4: Solve P M (ε) by using LINGO and obtain the solution and its corresponding objective vector. 5: end for 6: Remove dominated points and output the final solution set.

VII. COMPUTATIONAL EXPERIMENTS

In this section, benchmark instances in [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF] and largersize randomly generated instances are tested to evaluate the performance the proposed model algorithms. Both models P 0 and P are solved by using commercial optimization software LINGO (V14.0) and the proposed algorithms are coded in C++ with Visual Studio 2010. All the computational experiments are conducted on a PC with 2.5 GHz and 2.95 GB RAM under windows 10. The computational time (CPU time) spent by each method is limited to 10 min (600000 ms) for facilitating fast decision-making in actual emergencies.

A. Model Comparision

To evaluate the performance of the proposed model, we test the benchmark instances in [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF] by directly using P and P 0 solved by LINGO. For the sake of convenience, let Var and Cons denote the number of variables and constraints, respectively. Table II presents their comparison results.

From Table II, we can see that our model P employs much fewer variables and constraints than the existing one P 0 . P can exactly solve four of the five instances, while P 0 fails to solve any of them within 10 minutes. This indicates that our model is much more efficient than the existing one, because the proposed model is more compact and tighter than the one in [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF].

B. Real-Life Instance

We then test the real-life instance from the emergency scheduling of forest fires in Huzhong region located in Mt. Daxing'anling of Heilongjiang Province, China [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. In this instance, seven fire points are involved and 40 rescue vehicles are available for fighting fires. For more details, interested readers can refer to [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. Table III (called TRZ for short) and the proposed algorithms, in which columns 3 to 5 list the Pareto points. For each point, the first value represents the fire-extinguishing time, while the second one denotes the number of dispatched fire engines.

It can be seen in Table III that each method obtains a set of Pareto solutions. To better understand different possible decisions, take the first and last points of TRZ as an example, the results mean that it will take 6.08 h to all fire points with 40 fire engines for the former, while 39.6 h will be consumed with 29 fire engines for the latter. Decisionmakers can select a preferred solution from the Pareto front. Furthermore, from Table III we can conclude that:

1) In terms of solution quality, the number of Pareto-optimal solutions obtained by our algorithms is more than that by TRZ. Moreover, three new non-dominated points (i.e., (6.06, 40), (6.97, 38) and (15.45, 32)) are found by our algorithms. The heuristic can also obtain the Pareto front as the DP algorithm. This shows that the proposed heuristic is able to obtain high-quality solutions. 2) In terms of computational time, both the DP and heuristic algorithms are much more efficient in solving this reallife instance than TRZ and LINGO. Moreover, it can be observed that the heuristic performs better than the DP algorithm in terms of computational efficiency.

3) The fire-extinguishing time decreases as the number of dispatched fire engines increases due to the strictly decreasing relationship between them, which matches the conclusion made in Proposition 2.

C. Benchmark Instances

To further show the performance of the proposed methods, we conduct numerical comparison tests on a set of benchmark instances in [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF] with a different number of fire points, i.e., 10, 15, 20, and 50. Let R TRZ , R H ,a n dR DP denote the coverage rate of non-dominated points found by TRZ, the proposed DP and heuristic algorithms to the exact Pareto front, respectively. This indicator measures the solution quality. As the DP algorithm can obtain all the non-dominated points over the Pareto front, R DP is 100%. Note that both R TRZ and R H are equal to or less than 100%, as their obtained solutions are not guaranteed to be Pareto-optimal. For brevity, let |F | denote the number of obtained non-dominated points in the Pareto front. With the above notations, the comparison results are summarized in Table IV.

It can be seen from Table IV that the values of R H and R DP are 100% for all the instances, while that of R TRZ varies from It is worthwhile to point out that the value of R TRZ is 0 for the largest-size instance with n = 50 and m = 250, which implies that TRZ fails to derive Pareto-optimal solutions for this instance; whereas the value of R H does not decrease with the problem size. The results show that both the proposed DP and heuristic methods outperform TRZ in terms of solution quality. Furthermore, it can be seen that our algorithms can obtain the Pareto front for these instances in short time, i.e., the computational time spent by the proposed DP and heuristic algorithm slightly vary from less than 1 to 6 ms and from 12 to 82 ms, respectively, while the computational time of TRZ and LINGO varies from 525 to 51526 ms and 22000 to over 600000 ms, respectively, which sharply increases with the problem size. This shows that our algorithms are more efficient than TRZ and LINGO. In addition, the heuristic method achieves the best performance in terms of computational time.

In practical emergency scheduling, it is necessary to take into account the impact of the quantity of available rescue resources, i.e., the number of available fire engines to our problem. We also test the benchmark instances under different settings of the number of available fire engines (i.e., m). As it is very time-consuming to solve the model with LINGO, we only compare our algorithms with TRZ. Their comparison results are reported in Table V.

From Table V, we can see that the proposed DP and heuristic can obtain all non-dominated points (i.e., Pareto front), while TRZ can at most find 75% non-dominated points in the Pareto front. In addition, for each given n,t h ev a l u eo fR TRZ decreases with the value of m, while the coverage rates of our algorithms keep the same. The results show that the proposed algorithms are more effective than TRZ in solving these instances. Moreover, it can be seen that for a fixed n, the number of non-dominated points |F | increases with m, which implies the complexity of the problem increases with m. In terms of computational time, we can see that the proposed algorithms consume far less time than TRZ and also, the heuristic algorithm is the fastest one. The results again confirm

D. Larger-Size Instances

To further evaluate the performance of the proposed methods for solving larger-size problems, we perform numerical tests with more fires points varying from 100 to 1000. Note that although these large-size instances might not happen in real life, we conduct such experiments to further demonstrate the performance of the proposed algorithms. For each problem size, five instances are tested. All instances are randomly generated based on the way in [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF], which is described as follows. The distance among points in the network d ij is uniformly and randomly generated in the interval [START_REF] Han | Cuckoo search and particle filter-based inversing approach to estimating defects via magnetic flux leakage signals[END_REF]100]. The driving speed of the rescue vehicle from vertices i and j, ∀i, j ∈ V ′ , v ij is set to be 54 km/h. The fire-extinguishing speed of a fire engine v f is set as 2.5 m/min. The fire spread speed v si is uniformly and randomly generated in the interval [START_REF] Chiu | Real-time mobilization decisions for multipriority emergency response resources and evacuation groups: Model formulation and solution[END_REF][START_REF] Knott | Vehicle scheduling for emergency relief management: A knowledge-based approach[END_REF] m/min. As the TRZ consumes too much computational time for largesize problems, we only compare our two algorithms in this subsection. For the sake of convenience, let T DP and T H denote the computational time by the proposed DP and heuristic algorithms, respectively. The computational results are summarized in Table VI. Note that each value in the table is the average value of five instances.

From Table V, we can conclude that: 1) The value of R H is 100% for all problem sets, which implies that the proposed heuristic obtains the Pareto fronts of all instances.

2) The proposed algorithms obtain the Pareto fronts with up to 1000 fire points and 5000 fire engines within 79938 ms (less than 80 s). The results show that our algorithms are effective and efficient in solving the large-size instances. Also, it can be observed that T H is less than T DP over all problem sets, and T DP rapidly increases with the problem size, whereas T H increases in a much slower pace than DP does. This indicates that the heuristic method is more efficient than the DP algorithm. Besides, the results in this table suggest that |F | has an obvious increasing trend as the number of available fire engines (i.e., m) increases. This can be explained in part by Theorem 1.

E. Results of Multi-Depot Fire Engine Scheduling

We also test the real-life instance by considering multiple depots to evaluate the proposed model P M . Three cases with |K |=2, 3, and 4 are considered, where |K | denotes the number of fire service points (i.e., depots). According to [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF], the distances between the newly added depots and fire points are uniformly and randomly generated in [START_REF] Han | Cuckoo search and particle filter-based inversing approach to estimating defects via magnetic flux leakage signals[END_REF]100]. The total number of fire engines for each case is set the same as the work [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. From Table VII, we can see that the proposed algorithm (i.e., Algorithm 2) is able to obtain the Pareto front of the real-life instance with multiple depots within the given time. Moreover, the computational time increases as the number of depots increases, which is mainly because the number of variables and constraints increases with |K |. This indicates that the complexity of the problem increases with the number of depots. Besides, it can be observed that the number of non-dominated points found for each case is the same, which means that Algorithm 2 can find a non-dominated point at each iteration for each case. This implies that a non-dominated solution exists for each value of f 2 .

From Fig. 2, we can see a clear trade-off between the total rescue time and the number of fire engines dispatched for each |K |. The Pareto fronts derived can provide useful information for decision-makers to achieve a balance between rescue time and dispatched resources. In particular, decisionmakers can select their most preferred solution from the solution set. For instance, if they prefer a solution with the least emergency response time, then the last solution can be chosen.

VIII. CONCLUSIONS

This work has studied a bi-objective emergency scheduling problem for forest fires to simultaneously minimize the total fire-extinguishing time and the total dispatched fire engines. An improved bi-objective integer program containing less variables and constraints has been proposed for the studied problem. Experimental results indicate that the improved model significantly outperforms the existing one [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]. Based on the property analysis results of the problem, an exact dynamic programming algorithm and a fast greedy heuristic method are devised. Computational results for a real-life instance and benchmark instances show that the proposed DP and heuristic algorithms significantly outperform the existing method [START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF] in terms of both solution quality and computational efficiency. The Pareto fronts of all the benchmark instances are obtained. Larger-size randomly generated instances are also tested to evaluate the performance of the proposed algorithms. Experimental results indicate that large-size problems with up to to 1000 fire points and 5000 fire engines can be exactly solved within 80s on a personal computer. This demonstrates the effectiveness and efficiency of the proposed algorithms in solving large-scale rescue vehicle scheduling problems to deal with forest fires. Moreover, this study has also constructed a multi-depot fire engine scheduling model for the multi-depot case and the ε-constraint method is adapted to solve it. The trade-off between the two objectives to be optimized is clearly illustrated.

Although the proposed model is much more efficient than the existing one, it could not address the case with heterogeneous fire engines. One future research direction is to construct new and effective models for the heterogeneous fire engine scheduling. Other future research following this work may include extending the studied problem to more general cases. For example, in practice, the travel time from the fire service station to fire points can dynamically change with the roadway conditions. Extending the proposed model to a stochastic or time-dependent one by considering stochastic or varying travel time would be an important direction. Besides, it can be found that the multi-depot fire engine scheduling problem is much more difficult to solve. New effective and efficient approaches [START_REF] Tarantilis | A hybrid metaheuristic algorithm for the integrated vehicle routing and threedimensional container-loading problem[END_REF]- [START_REF] Li | Colored traveling salesman problem[END_REF] should be designed for addressing such optimization problems.

Fig. 1 .

 1 Fig. 1. Schematic illustration of the basic workflow of the forest fire-fighting.

  For |K |= 2 ,3 ,a n d4 ,w es e tm 1 = m 2 = 20; m 1 = 15, m 2 = 13, m 3 = 12; and m 1 = m 2 = m 3 = m 4 = 10, respectively. The other related data is kept unchanged. Table VII and Fig. 2 present the computational results.

TABLE I NOTATIONS

 I USED TO DESCRIBE THE PROBLEM

TABLE II MODEL

 II COMPARISON RESULTS FOR SOLVING INSTANCES IN[START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]

	TABLE III
	COMPARISON RESULTS FOR A REAL-LIFE INSTANCE

  reports the comparison results between the method proposed by Tian et al.[START_REF] Tian | Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm[END_REF]

TABLE IV COMPARISON

 IV RESULTS FOR BENCHMARK INSTANCES

TABLE V COMPARISON

 V RESULTS FOR BENCHMARK INSTANCESWITH DIFFERENT m VALUES 0 to 58.33%. This implies that proposed greedy heuristic can also obtain the Pareto front just like the DP algorithm, while TRZ fails to achieve it. Moreover, we can see that the value of R TRZ decreases as the problem size increases.

TABLE VI COMPUTATIONAL

 VI RESULTS FOR LARGER-SIZE RANDOM INSTANCES the efficiency of the proposed algorithms.

TABLE VII COMPUTATIONAL

 VII RESULTS FOR THE REAL-LIFE INSTANCE WITH MULTIPLE DEPOTS Fig. 2. The Pareto fronts of the real-case instance with multiple depots.

ACKNOWLEDGEMENT

We would like to thank the associate editor and four anonymous reviewers for their constructive comments and suggestions, which greatly improved the quality of this paper. We are also grateful to Dr. Guangdong Tian for providing his test instances, enabling us to compare the algorithms.