3D Multi-Drone-Cell Trajectory Design for Efficient IoT Data Collection

05/29/2019
by   Weisen Shi, et al.
Nanjing University
Xidian University
0

Drone cell (DC) is an emerging technique to offer flexible and cost-effective wireless connections to collect Internet-of-things (IoT) data in uncovered areas of terrestrial networks. The flying trajectory of DC significantly impacts the data collection performance. However, designing the trajectory is a challenging issue due to the complicated 3D mobility of DC, unique DC-to-ground (D2G) channel features, limited DC-to-BS (D2B) backhaul link quality, etc. In this paper, we propose a 3D DC trajectory design for the DC-assisted IoT data collection where multiple DCs periodically fly over IoT devices and relay the IoT data to the base stations (BSs). The trajectory design is formulated as a mixed integer non-linear programming (MINLP) problem to minimize the average user-to-DC (U2D) pathloss, considering the state-of-the-art practical D2G channel model. We decouple the MINLP problem into multiple quasi-convex or integer linear programming (ILP) sub-problems, which optimizes the user association, user scheduling, horizontal trajectories and DC flying altitudes of DCs, respectively. Then, a 3D multi-DC trajectory design algorithm is developed to solve the MINLP problem, in which the sub-problems are optimized iteratively through the block coordinate descent (BCD) method. Compared with the static DC deployment, the proposed trajectory design can lower the average U2D pathloss by 10-15 dB, and reduce the standard deviation of U2D pathloss by 56

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

05/30/2019

Multi-Drone 3D Trajectory Planning and Scheduling in Drone Assisted Radio Access Networks

Drone base station (DBS) is a promising technique to extend wireless con...
03/01/2021

Technical Report for A Joint User Scheduling and Trajectory Planning Data Collection Strategy for the UAV-assisted WSN

Unmanned aerial vehicles (UAVs) are usually dispatched as mobile sinks t...
08/18/2018

On Energy Efficiency of Networks for Composable Datacentre Infrastructures

This paper evaluates the optimal scale of datacentre (DC) resource disag...
07/23/2021

Trajectory Design for UAV-Based Internet-of-Things Data Collection: A Deep Reinforcement Learning Approach

In this paper, we investigate an unmanned aerial vehicle (UAV)-assisted ...
11/28/2019

Joint User Association and Resource Allocation in the Uplink of Heterogeneous Networks

This letter considers the problem of joint user association (UA), sub-ch...
03/20/2018

AC/DC: In-Database Learning Thunderstruck

We report on the design and implementation of the AC/DC gradient descent...
10/02/2019

A Machine Learning framework for Sleeping Cell Detection in a Smart-city IoT Telecommunications Infrastructure

The smooth operation of largely deployed Internet of Things IoT applicat...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

In the future Internet-of-things (IoT), the IoT devices are foreseen to be widely deployed to collect the data required by a myriad of applications. However, in certain areas where the terrestrial network coverage is unavailable due to signal blockage, spectrum scarcity or low IoT transmit power, collecting the IoT data becomes a challenging issue [1]. Although such an issue can be alleviated through densely deploying massive BSs or small-cells, the prohibitive capital expenditure (CapEx) and operating expenditure (OpEx) are unacceptable for IoT operators. In addition, the IoT data traffic is spatio-temporally unbalanced and dynamic-distributed in dedicated areas [2], which leads to frequent network congestion in the case of fix-deployed networks. To address the connection and flexibility shortages, the state-of-the-art drone communication technology is leveraged to support IoT data collections [3].

Various field experiments have verified the potential of drone, i.e., unmanned aerial vehicle (UAV), to support communication services for terrestrial users. Equipped with certain communication modules, the flying drone can perform as the drone cell (DC) which collects the data from IoT devices through the user (IoT device)-to-drone (U2D) wireless links. Compared with the terrestrial BS, DC has distinct advantages:

  • Line-of-sight (LoS) connection: User-to-BS (U2B) wireless links in terrestrial IoT are frequently blocked by the terrain or buildings, and such non-line-of-sight (NLoS) wireless links seriously degrade the U2B communications [4]

    . In contrast, the flying DCs are capable of adjusting 3D positions, which can ensure higher probability to connect terrestrial IoT devices via the highly reliable LoS links

    [3].

  • Dynamic deployment: Compared with terrestrial BSs built on fixed locations, DCs can change their deployments according to the spatio-temporal traffic variations, and assigned to different controllers or IoT devices on demands [5].

  • Fully-controlled mobility: The flying traces and behaviors of any DC are fully controlled by the network providers, which empowers DCs with the dynamic deployment feature and facilitates inter-DC collaborations [6].

Leveraging the advantages of DC, how to control the 3D placements of DCs to improve the network performance is of great importance yet very challenging due to the complicated 3D mobility of DC and dynamic DC-to-ground (D2G) link quality. In the literature, there have been several studies on the 3D placements of single/multiple DCs to support terrestrial IoT, which can be classified as two categories, i.e., static DC deployment and DC trajectory design. The static DC deployment investigates the optimal hovering positions of DCs to maximize the communication performance of associated users. Bor-Yaliniz

et al. proposed a bisection search based algorithm to deploy single DC [7]. Kalantari et al. optimized the multiple DCs static deployment through the swarm intelligence based algorithm [8] [9]. To maximize the IoT information collection gain, Mozaffari et al. leveraged the clustering based method to design optimal hovering positions for multiple DCs [10]. However, the static deployment can hardly guarantee the fairness among users. Particularly, the IoT devices located at the edge of the DC’s radio coverage can suffer severer pathloss compared with devices located at the center of the radio coverage, which leads to unbalanced communication performance.

Fig. 1: DC enabled IoT data collection.

To promote the fairness among IoT devices, DC trajectory design schemes are proposed in some recent works, in which DCs periodically fly over associated users via optimized trajectories, and communicate with dedicated users in scheduled slots [12]-[16]. Following pre-determined DC trajectories, Li et al. investigated the packets transmission relayed by multiple DCs from the sensor to BS [11]. Zeng et al. proposed a general framework for trajectory and communication optimization of a energy-efficient point-to-point UAV-ground communication system [12]. In [13], Wu et al.. formulated the non-convex DC trajectory design problem in which the user association, DC horizontal trajectory design and DC transmitting power control are jointly optimized; the delay constraint is further considered in the extension work focusing on delay-sensitive services [14]. However, there are some strong assumptions and model simplifications in current DC trajectory design studies. For instance, in current works, the flying altitudes of all DCs are usually set to be a constant value to simplify the optimization process (i.e., only 2D trajectory is optimized), and the adopted D2G channel model is the classic Friis model which fails to characterize the distinct features of DCs [13] [14]. Besides, most existing trajectory design works idealize or ignore the drone DC-to-BS (D2B) link quality constraints, which however is essential for drone-assisted networks.

To address the aforementioned issues, in this paper, we study the multiple DCs 3D trajectory design to facilitate the IoT data collection. We leverage the state-of-the-art U2D [15] and D2B [16] pathloss models, and the DC altitudes are taken as optimization variables to further improve the network performance. We formulate the 3D multi-DC trajectory design as a mixed integer non-linear programming (MINLP) problem, in which the summation of average U2D pathloss is minimized under the constraints of D2B link qualities, user fairness and horizontal/vertical speed limitations. Given the constant transmit power, bandwidth and interference-level of IoT devices, the lower average U2D pathloss indicates higher network throughput. As the formulated MINLP problem is very difficult to be solved directly, we decouple the original MINLP problem into multiple quasi-convex or linear sub-problems. Each sub-problem optimizes the user association, user scheduling, DC horizontal trajectories and flying altitudes, respectively. By adopting the block coordinate descent (BCD) method, we then propose a 3D multi-DC trajectory design algorithm in which those sub-problems are iteratively optimized to achieve the optimal result of the MINLP problem. According to the simulation results, our algorithm can reduce the average D2U pathloss by , and lower the standard deviation by 56% when compared with the static DC deployment scheme.

The remainder of this paper is organized as follows. In Section II, we introduce the system model. Then, we formulate the trajectory design problem for multiple DCs in Section III. In Section IV, we propose the 3D multi-DC trajectory design algorithm to solve the problem. Numerical results are shown in Section V, and the conclusion are given in Section VI.

Ii System Model

The overview of DC enabled IoT data collection is shown in Fig. 1, where multiple DCs are released to relay data uploading from IoT devices to one BS. We define the DC users as the IoT devices whose reliable U2B links are compromised due to blockage or resource shortage. The BS periodically detects the channel-state-information (CSI) of each U2B pair and assign compromised IoT devices into the DC user set . Given and the available DCs set , the BS implements the trajectory design algorithm to determine the optimal trajectory for each DC. DC users are considered as identical devices which transmit data to the BS with identical data collection period, bandwidth, and transmit power. Since the proposed trajectory design algorithm can schedule the transmit time of all users associated to one DC to prevent transmit-interval overlapping, and we assume each DC can distinguish its associated users’ signal efficiently from others, the inter-user interference issue is prevented in our system. The cardinalities and represent the number of DC users and available DCs, respectively.

Ii-a U2D and D2B Channel Models

Both U2D and D2B wireless links are modelled based on the recent D2G channel research [15] [16]. According to [15], we can express the U2D pathloss as

(1)

where is the U2D horizontal distance, and is the DC flying altitude. and are carrier frequency and speed of light in Hz and m/s, respectively. and are environment-dependent LoS and NLoS pathloss offsets, respectively. is the U2D LoS probability, which is expressed as

(2)

where and are environment-dependent parameters. To prevent the interference to U2B communications, as well as to provide additional spectrum resource for DC users, is expected to be different from the traditional IoT band, such as Wi-Fi or NB-IoT [17] bands. The average D2B pathloss is calculated by [16]

(3)

where and denote the D2B horizontal distance and the DC to BS antenna vertical angle, respectively. , , , , and are all environment-dependent parameters. Since the 850 MHz band is used for LTE transmissions [16], (3) contains no parameter representing carrier frequency.

Ii-B DC Trajectory Model

For one DC , we consider that it serves the associated user set through a time division multiple access (TDMA) mode. Within one trajectory period , DC flies over all its associated users and sequentially severs them according to the scheduling result. is set to the same value as the data collection period of users, in order that each user can transmit the collected data once within one period. By discretizing into equal-time slots, trajectory of DC within one period can be modeled as an

-length sequence composed by three-dimensional vectors:

(4)

where , and denote the 3D coordinates of DC at slot . is the set of slots in one period. In this work, we define when user is associated to DC for , otherwise

. The U2D communication scheduling is denoted by the binary variable

for . If user is severed by DC in slot , is set as ; otherwise, .

Several trajectory constraints are considered in our work: 1) The trajectory of each DC has to be a 3D closed curve since the DC must return to the start point by the end of each period for a new period. 2) Since the power consumption models for horizontal and vertical moving of DC are different, it is reasonable to set the maximum horizontal speed and maximal vertical speed , respectively. 3) In any slot , one DC can serve at most one user ; in all slots, one user can only be associated to one DC. 4) The slots amount scheduled to one within one period cannot be smaller than a pre-defined threshold , which indicates the minimal time to complete each U2B data transmission.

Iii Problem Formulation

We define as the coordinate of , as the 2D projection of DC on X-Y plane, and as the horizontal distance between and . Define the elevation angle between DC and BS antenna as in degree. is the Euclidean norm.

By substituting and into (1) and (3) respectively, we can calculate the U2D pathloss between DC and user in slot by

(5)

as well as the D2B pathloss between the BS and DC by

(6)

Assuming all users maintain their transmission power in each period, the effectiveness of the data collection for each U2D link is determined by the U2D pathloss. Therefore, we aim to minimize the total pathloss suffered by each DC in the trajectory design. Define the matrices , and . The 3D multi-DC trajectory design problem can be formulated as

(7a)
(7b)
(7c)
(7d)
(7e)
(7f)
(7g)
(7h)
(7i)
(7j)
(7k)

In (7), is the maximum number of users for one , and is the D2B pathloss threshold. (7c)-(7g) correspond to user association constraint 3) and 4) in section II. (7h)-(7j) represent the DC trajectory constraints 1) and 2) in section II. (7k) is the D2B pathloss constraint.

Due to the quadratic and exponential terms in (7a) and (7k), and the binary variable and , the 3D multi-DC trajectory design problem turns to be an MINLP problem [7]. Besides, the objective function (7a) and constraint (7n) are both non-convex for DC trajectory , which increases the difficulty to solve the problem.

Iv 3D Multi-DC Trajectory Design Algorithm

The non-convex problem (7) can be transformed into solvable forms (e.g. quasi-convex or LP) by assuming parts of the decision variables as constants. Then, we can decouple the MINLP problem into multiple sub-problems that are solvable for parts of the decision variables, and iteratively solve them by applying the BCD method. Specifically, we divide the decision variable in problem (7) into four blocks (i.e. , , and ), and formulate following sub-problems.

Given the pre-defined trajectories of multiple DCs (constant , and ), the user association sub-problem can be written as an integer linear programming (ILP) problem:

(8)

Problem (8) can be efficiently solved through the branch and bound method, which is well supported by Gurobi solvers.

Based on the optimized and pre-defined trajectory, the U2D communication scheduling sub-problem is also transformed to an ILP problem:

(9)

Same as problem (8), the Gurobi optimizer can efficiently solve problem (9).

According to (6), for given , and , (7a) is still non-convex for . However, by further decoupling into element optimization variable , (7a) is both quasi-convex and non-decreasing function to U2D horizontal distance , when other keep constants. Since (7a) can reach its minimal value when for any within its available range, the objective of minimizing (7a) equals minimizing , which is a quadratic convex function for . On the other hand, [18] proves that the available range of for constraint (7k) forms a convex set in X-Y plane. Therefore, the sub-problem to find optimal can be formulated as

(10)

whose optimal can be calculated by solving the following quadratic convex optimization problem:

(11)

For any given , and , the sub-problem optimizing is also non-convex to due to the -form and a -form terms. Same as preceding decoupling process of , can also be decoupled into element optimization variable , and the sub-problem to optimize is

(12)

Problem (12) can only be proved quasi-convex to when other keep constant. However, since is the summation of one non-increasing function of and one non-decreasing function of , we can argue that problem (12)’s objective function has only one global minimum. To achieve the optimal which minimizes problem (12), we can leverage the Newton-Raphson method by iteratively calculating the following function:

(13)

where the iteration stops when and .

The sub-problems of problem (7) can all be optimized respectively with other optimization variables keeping constant. Therefore, problem (7) can be solved through iteratively optimizing those sub-problems until the results converge, which yields to the classic BCD method. Based on the BCD method, we propose the 3D multi-DC trajectory algorithm as shown in Algorithm 1. , , denote the multiple DCs’ user association, U2D communication scheduling and DC trajectories calculated after iteration , respectively. Since the global optimal results of each sub-problem can be achieved accurately, the proposed algorithm ensures convergence [13] [19].

1:Initiate DC set and their initial trajectory composed by and .
2:Initiate initial U2D communication scheduling .
3:, = .
4:while  do
5:     Solve (8) to obtain with , and .
6:     Solve (9) to obtain with , and .
7:     for  do
8:         Solve (10) for with , and .
9:         Update with .
10:     end for
11:     for  do
12:         Solve (12) for with , and .
13:         Update with .
14:     end for
15:     Update with and .
16:     .
17:     .
18:end while
Algorithm 1 3D multi-DC trajectory design algorithm
(a) 3D view
(b) Vertical view
(c) Altitudes
Fig. 2: Trajectory design results of 5 DCs to cover 20 users.

V Simulation Results

We implement the proposed algorithm in simulations by using Gurobi solver [20]. The parameters of both U2D and D2B pathloss models are set as suburban scenario [16]. The carrier frequency for U2D communication is set as , which is widely supported by commercial drone products and prevents interference to the cellular system. D2B communications use the LTE band [16]. As the initial trajectories, DCs are uniformly deployed over the radio coverage area of terrestrial BS with the same flying altitude . Without loss of generality, we treat one slot as the minimal time unit to calculate related variables including , and . According to the general specifications of commercial drones whose maximal horizontal and ascent/descent speeds are and respectively, the approximate value of is around . Table. I summarizes the detail simulation parameters.

Parameter Name Value
BS radio coverage radius
User(IoT device) number
U2D parameters (0.1,21,4.88,0.43)
D2B parameters (3.04,-23.29,-3.61,4.14,20.7)
Carrier frequencies (U2D, D2B)
Duration of one period
D2B pathloss constraint
Minimal number of slots
Maximal horizontal speed
Maximal vertical speed
Trajectory difference for each slot
TABLE I: Simulation Parameters

Fig. 2 shows an example of trajectory design result in which five DCs relay data from 20 IoT users. The closed curves with different markers in Figs. 2(a) and 2(b) denote the trajectories of different DCs, while the squares on the X-Y plane represent the users. Users are associated to their corresponding DCs with same colors. Fig. 2(c) shows the variations of flying altitude within one period, where the lower bound of all trajectories near is constrained by the D2B link quality constraint. As shown in Figs. 2(a) and 2(b), for each DC, the optimized trajectory can fly over all its associated user and form a closed trajectory in 3D space. In Fig. 2, the time interval between two adjacent dots on the same trajectory is exactly one slot.

Fig. 3: Average U2D pathloss comparison between trajectory design and static deployment.

In Fig. 2(b), we can see that more than dots on each trajectory are overlapped above the associated users positions. By joint analyzing Figs. 2(b) and 2(c), we can justify that DCs are prone to hovering above the associated users, while spending less slots for the travel process between two hovering positions. Such a “hovering effect” in the optimized trajectories can indicates the effectiveness of our proposed algorithm in minimizing average U2D pathloss. Because with given flying altitude, the average U2D pathloss can be minimized when U2D horizontal distance is zero, i.e., hovering above the user.

Since most existing trajectory planning results are based on different D2G channel models with fixed-height (i.e. [13]

), we choose one static DC deployment algorithm, the per-drone iterated particle swarm optimization (DI-PSO) algorithm

[18], as the benchmark to highlight the efficiency of our proposed algorithm. The average U2D pathloss performance achieved by the proposed 3D multi-DC trajectory design algorithm, as well as the static DC deployment scheme are compared in Fig. 3 and Fig. 4. In both figures, the blue curves represent static DC deployment while the red curves represent DC trajectory design.

Fig. 4: U2D pathloss CDF curves of trajectory design and static deployment.

From Fig. 3, we can see that the U2D pathloss performance of both algorithms are improved as the available DC number increases. However, the average U2D pathloss of our trajectory design solution remains less than that of the DI-PSO solution by . The CDF curves in Fig. 4 further confirm that the trajectory design algorithm can maintain the U2D pathloss below given thresholds with higher probability when compared with static DC deployment.

The user fairness promotion provided by the DC trajectory design algorithm is indicated by the error bars in Fig. 4 and the U2D pathloss standard deviation comparison in Table II. and are U2D pathloss standard deviations for DC trajectory design and static DC deployment, respectively. In Fig. 4, we can observe that the average U2D pathloss of static deployment ranges from to , which is four times as much as that of the trajectory design algorithm. From Table II, we can conclude that the trajectory design can reduce the U2D pathloss standard deviation by on average compared with static deployment.

DC number
TABLE II: U2D pathloss standard deviation comparison

Vi Conclusion

In this paper, we have investigated the 3D multi-DC trajectory design for efficient IoT data collection. An MINLP problem has been formulated to minimize the summation of average U2D pathloss. To solve the MINLP problem, we have decoupled it and formed multiple sub-problems in which the user association, U2D communication scheduling, horizontal trajectories, and flying altitudes are optimized, respectively. Leveraging the BCD method, we have devised the 3D multi-DC trajectory design algorithm to solve the MINLP problem by solving sub-problems iteratively. Simulation results have shown that the proposed DC trajectory design algorithm can achieve average U2D pathloss reduction, and promote pathloss standard deviation by more than when compared with the static DC deployment. In future works, we will analyze the impacts of initial deployments, horizontal and vertical flying speeds, as well as inter-DC safe distance on the trajectory design, and investigate the communication and computation resources allocation of multiple DCs with the optimized trajectories.

Acknowledgment

This work is supported by the National Natural Science Foundation of China under Project 91638204 and the Natural Sciences and Engineering Research Council (NSERC), Canada.

References

  • [1] N. Cheng, W. Xu, W. Shi, Y. Zhou, N. Lu, H. Zhou, and X. Shen, “Air-ground integrated mobile edge networks: Architecture, challenges, and opportunities,” IEEE Commun. Mag., vol. 56, no. 8, pp. 26–32, 2018.
  • [2] M. Tao, K. Ota, and M. Dong, “Locating compromised data sources in IoT-enabled smart cities: A great-alternative-region-based approach,” IEEE Trans Ind. Informat., vol. 14, no. 6, pp. 2579–2587, 2018.
  • [3] M. Mozaffari, W. Saad, M. Bennis, Y.-H. Nam, and M. Debbah, “A tutorial on UAVs for wireless networks: applications, challenges, and open problems,” 2018. [Online]. Available: https://arxiv.org/abs/1803.00680.
  • [4] J. Liu, Y. Shi, Z. M. Fadlullah, and N. Kato, “Space-air-ground integrated network: A survey,” IEEE Commun. Surveys Tuts., early access, 2018.
  • [5] Z. Xue, J. Wang, G. Ding, H. Zhou, and Q. Wu, “Maximization of data dissemination in UAV-supported internet of things,” IEEE Wireless Commun. Lett., early access, 2018.
  • [6] S. Zhang, W. Quan, J. Li, W. Shi, P. Yang, and X. Shen, “Air-ground integrated vehicular network slicing with content pushing and caching,” IEEE J. Sel. Areas Commun., early access, 2018.
  • [7] R. I. Bor-Yaliniz, A. El-Keyi, and H. Yanikomeroglu, “Efficient 3D placement of an aerial base station in next generation cellular networks,” in Proc. IEEE ICC, Kuala Lumpur, Malaysia, May 2016, pp. 1–5.
  • [8] E. Kalantari, H. Yanikomeroglu, and A. Yongacoglu, “On the number and 3D placement of drone base stations in wireless cellular networks,” in Proc. IEEE VTC-Fall, Montreal, QC, Sept 2016, pp. 1–6.
  • [9] E. Kalantari, M. Z. Shakir, H. Yanikomeroglu, and A. Yongacoglu, “Backhaul-aware robust 3D drone placement in 5g+ wireless networks,” in Proc. IEEE ICC Workshops, Paris, France, May 2017, pp. 109–114.
  • [10] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Mobile unmanned aerial vehicles (UAVs) for energy-efficient internet of things communications,” IEEE Trans. Wireless Commun., vol. 16, no. 11, pp. 7574–7589, 2017.
  • [11] K. Li, W. Ni, X. Wang, R. P. Liu, S. S. Kanhere, and S. Jha, “Energy-efficient cooperative relaying for unmanned aerial vehicles,” IEEE Trans. Mobile Comput., vol. 15, no. 6, pp. 1377–1386, 2016.
  • [12] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” IEEE Trans. Wireless Commun., vol. 16, no. 6, pp. 3747–3760, 2017.
  • [13] Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for multi-UAV enabled wireless networks,” IEEE Trans. Wireless Commun., vol. 17, no. 3, pp. 2109–2121, 2018.
  • [14] Q. Wu and R. Zhang, “Common throughput maximization in UAV-enabled ofdma systems with delay consideration,” 2018. [Online]. Available: https://arxiv.org/abs/1801.00444
  • [15] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” IEEE Wireless Commun. Lett., vol. 3, no. 6, pp. 569–572, 2014.
  • [16] A. Al-Hourani and K. Gomez, “Modeling cellular-to-UAV path-loss for suburban environments,” IEEE Wireless Commun. Lett., vol. 7, no. 1, pp. 82–85, 2018.
  • [17] M. E. Soussi, P. Zand, F. Pasveer, and G. Dolmans, “Evaluating the performance of eMTC and NB-IoT for smart city applications,” in Proc. IEEE ICC, Kansas City, MO, May 2018, pp. 1–7.
  • [18] W. Shi, J. Li, W. Xu, H. Zhou, N. Zhang, S. Zhang, and X. Shen, “Multiple drone-cell deployment analyses and optimization in drone assisted radio access networks,” IEEE Access, vol. 6, pp. 12 518–12 529, 2018.
  • [19] D. P. Bertsekas, Nonlinear programming.   Belmont, MA, USA: Athena scientific, 1999.
  • [20] “Gurobi optimizer 8.1.” [Online]. Available: http://www.gurobi.com