Research Article | Open Access
Multiple Target Localization with Bistatic Radar Using Heuristic Computational Intelligence Techniques
We assume Bistatic Phase Multiple Input Multiple Output radar having passive Centrosymmetric Cross Shape Sensor Array (CSCA) on its receiver. Let the transmitter of this Bistatic radar send coherent signals using a subarray that gives a fairly wide beam with a large solid angle so as to cover up any potential relevant target in the near field. We developed Heuristic Computational Intelligence (HCI) based techniques to jointly estimate the range, amplitude, and elevation and azimuth angles of these multiple targets impinging on the CSCA. In this connection, first the global search optimizers, that is,are developed separately Particle Swarm Optimization (PSO) and Differential Evolution (DE) are developed separately, and, to enhance the performances further, both of them are hybridized with a local search optimizer called Active Set Algorithm (ASA). Initially, the performance of PSO, DE, PSO hybridized with ASA, and DE hybridized with ASA are compared with each other and then with some traditional techniques available in literature using root mean square error (RMSE) as figure of merit.
It is easy to deal with the far field targets because all the waves emitted by them are assumed to be plane waves and, by incorporating such assumption, one can characterize the signal model only by the Direction of Arrival (DOA) of the targets [1, 2]. On the other hand, it is comparatively difficult and challenging to localize the targets located close to the sensors array, that is, in near field. In this case, the plane wave assumption is no longer applicable and one has to deal with spherical waves. Thus, to localize the near field targets, the range information must also be incorporated into the signal model along with DOA [3, 4]. Near field targets localization, especially 3D (range, elevation angle, and azimuth angle), plays significant role in radar, cognitive radio networks, and array signal processing, since it is a preliminary step for adaptive beamformer to guide the main beam in preferred direction and simultaneously manage the nulls in the direction of jammers [5–7].
Several algorithms are available to address the issue of 3D near field source localization problem. In [8–10] different schemes are discussed but they not only are computationally expensive but also have the problem of pair matching between elevation and azimuth angles. To solve the pair matching problem, a two-stage separated steering vector-based algorithm is proposed in  but it has higher MSE and is computationally expensive as it requires more than 400 snapshots to achieve the results. It also fails to estimate the amplitude of sources which is also sometimes an important parameter to be estimated. Clearly the goal is to develop a scheme which must be able to jointly estimate the amplitude, range, elevation angle, and azimuth angle, should provide an improved MSE, and finally should be free of pair matching problem.
In order to achieve the desired goals Heuristic Computational Intelligence (HCI) techniques are available. In the last two decades, these techniques were being widely applied to a mixture of problems ranging from handy applications in industry and commerce to leading scientific research [12–14]. These techniques include Differential Evolution (DE), Particle Swarm Optimization (PSO), and Genetic Algorithm (GA). Besides the other useful properties, one of the most significant aspects of these techniques is their property of hybridization with any other heuristic or nonheuristic techniques. Due to this property, their reliability and applicability increase even more. In , an efficient scheme based on GA hybridized with Pattern Search (PS) is developed for the joint estimation of amplitude and 2D DOA of far field sources impinging on L shape arrays (1-L and 2-L shape arrays). Similarly, another useful algorithm is designed in  which is based on GA hybridized with Interior Point Algorithm (IPA) to jointly estimate the amplitude, range, and elevation angle of near field sources impinging on Uniform Linear Array (ULA).
In this paper, we assume a scenario of Bistatic Phase Multiple Input Multiple Output (MIMO) radar with colocated antennas which is capable of sending multiple beams with diverse wave forms on diverse targets. But, before getting into the mode of tracking targets by sending high gain sharp beams on the targets, it is required to estimate their spherical coordinates, that is, range and elevation and azimuth angles . To achieve this, the transmitter sends a coherent signal beam using any subarray of the transmitter such that a beam has a fairly wide solid angle so as to cover up any potential target in the space. The receiver in the first mode of this scenario is a CSCA of passive sensors consisting of two and subarrays. The reflected signals from different targets will impinge on these - arrays. Since we are interested in the ranges of the targets apart from 2D DOA, we treat the targets as near field in order to localize them completely. To jointly estimate the amplitude, range, and elevation and azimuth angles of the targets, we employ(a)PSO,(b)DE,(c)PSO hybridized with Active Set Algorithm (ASA),(d)DE hybridized with ASA.
In both hybrid schemes, the best individual results achieved through DE and PSO are given as starting point to ASA for further tuning. The performance criterion is made on the basis of mean square error (MSE) which is used as fitness function derived from maximum likelihood principle (MLP) . This fitness function is able to automatically pair the amplitude, range, and elevation and azimuth angles of a specific target. The other properties of this fitness function are its ease in concept, ease in implementation, robustness against noise, and requirement of single snapshot. The performances of PSO, DE, Particle Swarm Optimization hybridized with Active Set Algorithm (PSO-ASA), and Differential Evolution hybridized with Active Set Algorithm (DE-ASA) are not only compared with each other but also with some traditional techniques available in literature [10, 11] using RMSE as figure of merit.
The rest of the paper is organized as follows: Section 2 evaluates the problem formulation, while the proposed HCI techniques are discussed in Section 3. In the same way, results and discussion are carried out in Section 4, while Section 5 is dedicated to conclusion and future work directions.
2. Data Model for Near Field Targets
In this section, we developed a data model for near field targets impinging on CSCA placed on the receiver of Bistatic radar. The CSCA is composed of two symmetric sub-ULAs placed along -axis and -axis, respectively, as shown in Figure 1. Each ULA carries passive sensors, while the reference sensor is common for them. For , the data model at th and th sensor in the -axis and -axis subarray, respectively, can be represented asIn (1) and (2), is the propagation delay for th target between th and th sensors, while is the propagation delay for th target between th and th sensors. For th near field target, the signal phases and are parameterized mathematically as intermediate parameters in the form given asIn (3), , , , and are the function of th source amplitude , elevation angle , azimuth angle , and range . Mathematically, they can be expressed as In (1) and (2), and are Additive White Gaussian Noise (AWGN) added at th and th sensors, respectively. In vector form, the signal model can be represented aswherewhereNow clearly the problem in hand is to accurately and jointly estimate the unknown parameters (amplitudes, ranges, and elevation and azimuth angles) of the reflected signals from targets.
3. Heuristic Computational Intelligence
Heuristic Computational Intelligence (HCI) is a subfield of artificial intelligence that can be used for continuous and combinatorial optimization problems. The entire HCI algorithms are global optimization based methods which have stochastic or meta-heuristic characteristics which can be used for black box or gradient free problems. Basically, these kinds of algorithms make use of the iterative progress in a population and, to achieve the desired goal, the population is selected using parallel processing in random guided manner . The major advantages of HCI over traditional optimization methods include wide range of applicability, conceptual simplicity, knowledge usage capability, hybridization with any other scheme, robustness against dynamic changes, and ability of self-optimization . In literature, several HCI based algorithms are available but we are focusing on DE and PSO.
3.1. Particle Swarm Optimization (PSO)
PSO algorithm is inspired from bird flocking and school fishing which was proposed by Kennedy and Eberhart in 1995 . It can handle multimodal, discontinuous, and nonconvex optimization problems because it is a pseudooptimization method [20, 21]. The schematic diagram of PSO is shown in Figure 2, while its execution steps are given below.
Step 1 (initialization). Initialize the swarm randomly, that is, randomly producing particles. In this problem, each particle consists of genes, where represents the total number of available targets. Mathematically, it can be written asThe upper and lower bounds of , , , and are given aswhere and represent the lower and upper limits of ranges, while and are the upper and lower limits of amplitudes.
Step 2 (fitness function). We have used MSE as a fitness function to compute the fitness of every particle. For th particle, it can be mathematically written aswhere is called MSE which is derived from maximum likelihood principle . The total MSE for th particle can be mathematically written as where and are the number of sensors placed along -axis and -axis, respectively. In (12) and (13), and are given in (1) and (2), respectively, while, and can be defined asNow, store each particle as local best while the one having maximum fitness function is stored as global best particle.
Step 3a (update particle velocity). We used the following relation to update the velocity of each particle:In the right hand side of (15), the first term is called inertia or momentum which shows the previous velocity. The second term shows the private thinking and is also called as cognitive term of the particle while the last term explains the collective behavior of the population and is known as a social component. At the beginning, the value of to ensure that more weightage is given to local intelligence at the beginning; then there is gradual increase in the value of towards 0.9 which means that more weightage is given to collective intelligence at the end. In (15), both and are positive constants and for the ongoing problem .
The velocity is doubly bounded; that is,;;. and are taken to be 2 and −2, respectively, in order to avoid the premature convergence.
Step 3b (update particle position). The position of each particle can be updated as
Step 4a (choose local best particle). If fitness , replace with .
Step 4b (choose global best particle). Similarly replace the previous global best particle with as a new global best particle if fitness .
Step 5 (termination). Terminate, if any of the following conditions is satisfied; otherwise go to Step a:(i)The predefined total number of iterations/flights is executed for the algorithm.(ii)The desired MSE is attained which is .
Step 6 (hybridization). In this step, the best particle achieved through PSO is given to ASA for further improvement.
Step 7 (storage). Store the results for later on discussion and comparison.
3.2. Differential Evolution (DE)
DE was introduced by Storn and Price in 1997  which iteratively searches large spaces of candidate solution that improve the candidate solution with respect to specified measurement of quality. In other words, DE deals with the optimization problem in the following way: first it achieves the candidate solution and then, by using simple formula, it generates new solution by combining with the existing ones. As a result only those candidate solutions will survive which has best fitness or score. The block diagram of DE is shown in Figure 3, while its steps are given as follows.
Step 1 (initialization). The first step is similar to the one given for PSO in (8).
Step 2 (updating). Update all particles of the current generation “ge” from 1 to . Let us choose th particle from (8), that is, , where and , while “ge” represents the particular generation. The main goal of DE is to find out the particles of next generation, that is, , by adopting the following steps:(a) Mutation: select any three different numbers (particles) from 1 to , that is, , by keeping the following conditions: where “” is a constant whose value usually lies in the range 0.5 to 1.(b) Crossover: the crossover can be performed as where and is between 1 and chosen at random.(c) Selection operation: in this step, select the particle of next generation by using the following criterion of selection: where the and are defined in (8). Repeat this for all particles.
Step 3 (stopping criteria). The stopping criteria of DE are based on any of the following conditions being fulfilled:(I), where is a very small positive number OR.(II)The number of iterations is completed.Else go back to Step .
Step 4 (hybridization). In order to further refine the results of DE, we give the best particle to ASA as starting point.
3.3. Active Set Algorithm (ASA)
Active set methods are recursive procedures for solving sequence of equality-constrained quadratic subproblems effectively from majority of algorithms based on interior methods and simplex techniques due to the fact that they operate in two phases; first they focus on feasibility, while others focus on optimality. The applications in which their dominance is well established are quadratic programming including portfolio analysis, structural analysis, and optimal control and optimal mesh refinement in ODE or PDE constrained problems, sparse linear programming problems, and box constrained optimization problem [23–25]. Active set quadratic programming is exploited in many state-of-the-art optimization solvers such as QPOPT , SQOPT , and QPA (part of the GALAHAD library) . Although the global optimization methods (PSO and DE) alone can produced excellent results, their rate of convergence usually decreases with the increase of iteration especially in nonlinear problems. Therefore, to get accurate results, we need an efficient local search optimizer like ASA. The ASA is mainly used in constrained optimization problems whose basic purpose is to transform the problem into an easier solvable problem . The standard working procedure of ASA in the form of pseudocode is given as follows, while its flow diagram in terms of processing blocks is shown in Figure 4.
Generic Pseudocode of Active Set Algorithm
Start. Initialize ASA algorithm with random assignment of parameters.Step 1.Determine a suitable starting point.Step 2.Repeat till an optimal condition is achieved.Step 3.Solve the equality problem defined roughly by ASA.Step 4.Calculate the Lagrange multipliers for ASA and remove all those constraints having negative Lagrange multipliers.End Repeat
The mathematical formulation of ASA algorithm is given as.
The constrained optimization problem can be mathematically written aswhere and represent the lower and upper bounds, respectively, of the desired vectors. In early methods, penalty functions were used to transform the constrained problem into unconstrained one which is beyond or near the constrained boundary. In such a way, the constrained problems were solved by a sequence of parameterized unconstrained optimizations, which is obviously an inefficient way, that are now replaced by methods which focus on the solution of Karush-Kuhn-Tucker (KKT) equations such as ASA. The KKT equations are the essential condition for constrained optimization problem. In order to achieve the global solution point, the KKT equations are the sufficient and essential conditions to be satisfied. The KKT equations can be mathematically written aswhere , , are called Lagrange multipliers. Equation (21) represents the gradient cancelation of active constraints and the objective function at the solution point where the Lagrange multipliers are necessary to balance the deviation in magnitude of the constrained gradients and objective function.
For ASA, we used a MATLAB built-in optimization toolbox where the parameter setting is provided in Table 1.
4. Results and Discussion
In this section, several simulations are carried out to validate the performance of the proposed techniques based on PSO, DE, PSO-ASA, and DE-ASAs. For theoretical consideration among the optimization techniques, the strength of memetic computing techniques based on PSO-ASA and DE-ASAs is superior to the rest. Comparison of the results of these nonconventional meta-heuristic techniques is presented here in number of graphical and numerical illustrations to establish memetic optimization methodologies in terms of estimation accuracy, reliability, robustness, and convergence.
Simulation studies are present here in two parts; in part 1, we compared the performances of PSO, DE, PSO-ASA, and DE-ASA with each other in terms of estimation accuracy and convergence rate for different number of targets. In part 2, the performances of the two best techniques among them are compared with existing traditional algorithms [10, 11] by using MSE as a figure of merit. Every time, the number of sensors in both subarrays is taken to be same. In order to avoid the phase ambiguity, the interelement spacing between the consecutive elements is required . All the signals reflected back from targets are assumed to be statistically independent and having constant frequency. The received data at the output of each sensor are polluted by zero mean, unit variance AWGN. The values of elevation and azimuth angles are taken in degrees, while the values of ranges are taken in terms of wavelength . All the simulations for PSO and DE are carried out for 1000 independent trials.
4.1. Estimation Accuracy
In this subsection, estimation accuracy of PSO, DE, PSO-ASA, and DE-ASA is discussed without adding any noise to the system. In Tables 2–4, the estimation accuracy of all schemes is provided for two, three, and four sources, respectively, for different values of elevation angle, azimuth angle, range, and amplitudes. For this, the CSCA consists of 9, 13, and 17 sensors for 2, 3, and 4 targets, respectively. One can see from Tables 2–4 that each scheme produced better estimation accuracy for two targets but their accuracy degraded for three and four targets due to the presence of local minima. From these results, the advantages of hybridization are quite obvious as the estimation accuracy of PSO alone is less as compared to DE alone but its accuracy increases than DE when it is hybridized with ASA. Overall the DE-ASA produced better estimation accuracy as compared to the other three techniques. The second and third best results are given by PSO-ASA and DE alone, respectively.
In Table 5, mean, variance, and standard deviation are carried out for all schemes in terms of RMSE. For this simulation, two targets are considered to be impinging on CSCA having nine sensors. It is apparent that again DE-ASA produced better results as compared to the other algorithms while the second and third best results are once again given by PSO-ASA and DE alone.
In this section, the convergence of each scheme for different number of targets is discussed. By convergence we mean the total number of times a particular scheme got its desire results. As shown in Figure 5, the convergence of each scheme decreases with the increase of unknowns (targets) in the problem. However, the convergence of hybrid schemes are less degraded and they maintained fairly good convergence every time. The first and second best schemes among them are DE-ASA and PSO-ASA as they maintained better convergence. The third best convergence is given by DE alone. From the above discussion, the first and second best schemes are DE-ASA and PSO-ASA, respectively, so, to summarize the discussion, we shall be limited to the discussion of these two hybrid schemes only in the upcoming discussion.
4.3. Proximity Effects
In this subsection, the results are carried out to assess the proximity effects of elevation and azimuth angles for DE-ASA and PSO-ASA. We considered 3 sources impinging on CSCA, composed of 13 sensors. As provided in Tables 6 and 7, both of the hybrid schemes produced fairly good estimation accuracy and convergence for closely spaced elevation and azimuth angles. However, the DE-ASA produced better results as compared to PSO-ASA.
4.4. Estimation Accuracy for DOA on Reference Axis
Some of the elevation angles and azimuth angles are considered to be critical angles where the performance of most algorithms degraded. We considered the three targets at , , and . In Figure 6(a), the estimation accuracy of both hybrid schemes for elevation angles is evaluated on reference axis ( and ). It can be seen that both hybrid schemes produced significant error ( to ) especially for . However, the DE-ASA produced comparatively less errors as compared to PSO-ASA.
In Figure 6(b), the estimation accuracy of both hybrid schemes is evaluated for azimuth angles at , , and . In this case, both schemes are up to the mark and produced negligible error. Once again the DE-ASA produced better results as compared to PSO-ASA.
4.5. Comparison with Other Techniques Using Root Mean Square Error (RMSE)
In this subsection, we compared the RMSE of DE-ASA and PSO-ASA with existing nonheuristic techniques [10, 11] in the presence of noise. We considered two sources, and , where the CSCA consist of 17 sensors. In Figures 7(a)–7(c), both hybrid schemes produced lower RMSE as compared to the algorithms described in [10, 11]. The RMSE of both hybrid schemes is lower for the target located near to the array as compared to the target located comparatively away from the array as shown in Figure 7(c). Among all of them, the DE-ASA produced lower RMSE while the second best RMSE is obtained by the PSO-ASA scheme.
One of the other features of both hybrid schemes is that they can also be used for the amplitude estimation which is missing in [10, 11]. In Figure 7(d), RMSE is shown against SNR for the amplitude of the targets. Once again the DE-ASA scheme maintained a lower RMSE for both targets.
5. Conclusion and Future Work Direction
In this work, we developed Heuristic Computational Intelligence techniques to estimate jointly the amplitude, range, and elevation and azimuth angle of multiple targets impinging on the passive CSCA. In these computational techniques, first Particle Swarm Optimization (PSO) and Differential Evolution (DE) were used alone and then, to improve the results further, they were hybridized with ASA. It has been found that the hybrid DE-ASA schemes produced better results as compared to PSO and DE in terms of estimation accuracy and convergence and have also shown better results as compared to nonheuristic techniques in terms of RMSE.
Besides all the advantages, the proposed schemes have also few limitations:(i)All the schemes fail if the number of targets increases to that of the number of sensors in the CSCA.(ii)The performance of the proposed schemes degraded when the number of sources is more than four.In the future, one can use such kind of techniques in the following areas:(i)This kind of model for finding range, amplitude, and elevation angle can effectively be used for spectrum sensing in the cognitive radio domain in which the source shall be mainly primary user station. With each snapshot, one can find these updated parameters.(ii)One can use such kind of algorithms for null steering, sidelobe reductions, and main beam steering.
Conflict of Interests
All the authors of this paper declare that there is no conflict of interests regarding its publication.
- S. O. Al-Jazzar, A. Muchkaev, A. Al-Nimrat, and M. Smadi, “Low complexity and high accuracy angle of arrival estimation using eigen value decomposition with extension to 2D AOA and power estimation,” EURASIP Journal on Wireless Communications and Networking, vol. 2011, article 123, 2011.
- R. Roy and T. Kailath, “ESPRIT—estimation of signal parameters via rotational invariance techniques,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 7, pp. 984–995, 1989.
- J. Liang, X. Zeng, B. Ji, J. Zhang, and F. Zhao, “A computationally efficient algorithm for joint range-DOA-frequency estimation of near-field sources,” Digital Signal Processing, vol. 19, no. 4, pp. 596–611, 2009.
- Y. Wu, L. Ma, C. Hou, G. Zhang, and J. Li, “Subspace-based method for joint range and DOA estimation of multiple near field sources,” Signal Processing, vol. 86, no. 8, pp. 2129–2133, 2006.
- S. Darzi, T. S. Kiong, M. T. Islam, M. Ismail, S. Kibria, and B. Salem, “Null steering of adaptive beamforming using linear constraint minimum variance assisted by particle swarm optimization, dynamic mutated artificial immune system, and gravitational search algorithm,” The Scientific World Journal, vol. 2014, Article ID 724639, 10 pages, 2014.
- T. S. Kiong, S. B. Salem, J. K. S. Paw, K. P. Sankar, and S. Darzi, “Minimum variance distortionless response beamformer with enhanced nulling level control via dynamic mutated artificial immune system,” The Scientific World Journal, vol. 2014, Article ID 164053, 9 pages, 2014.
- S. Leng, W. Ser, and C. C. Ko, “Adaptive beamformer derived from a constrained null steering design,” Signal Processing, vol. 90, no. 5, pp. 1530–1541, 2010.
- N. Kabaoglu, H. A. Cirpan, E. Cekli, and S. Paker, “Deterministic maximum likelihood approach for 3-D near field source localization,” International Journal of Electronics and Communications, vol. 57, no. 5, pp. 345–350, 2003.
- J. Liang, S. Yang, J. Zhang, L. Gao, and F. Zhao, “4D near-field source localization using cumulant,” EURASIP Journal on Advances in Signal Processing, vol. 2007, Article ID 017820, 10 pages, 2007.
- R. N. Challa and S. Shamsunder, “Passive near-field localization of multiple non-Gaussian sources in 3-D using cumulants,” Signal Processing, vol. 65, no. 1, pp. 39–53, 1998.
- J. Liang, D. Liu, X. Zeng, W. Wang, J. Zhang, and H. Chen, “Joint azimuth-elevation/(-range) estimation of mixed near-field and far-field sources using two-stage separated steering vector-based algorithm,” Progress in Electromagnetics Research, vol. 113, pp. 17–46, 2011.
- D. B. Fogel, Evolutionary Computation: Toward a New Philosophy of Machine Intelligence, IEEE, Piscataway, NJ, USA, 2nd edition, 2000.
- D. B. Fogel and L. J. Fogel, “Using evolutionary programming to schedule tasks on a suite of heterogeneous computers,” Computers and Operations Research, vol. 23, no. 6, pp. 527–534, 1996.
- P. J. Angeline, G. M. Saunders, and J. B. Pollack, “An evolutionary algorithm that construct recurrent neural networks,” IEEE Transactions on Neural Networks, vol. 5, no. 1, pp. 54–64, 1994.
- F. Zaman, I. M. Qureshi, A. Naveed, J. A. Khan, and R. M. A. Zahoor, “Amplitude and directional of arrival estimation: comparison between different techniques,” Progress in Electromagnetics Research B, vol. 39, pp. 319–335, 2012.
- F. Zaman, I. M. Qureshi, A. Naveed, and Z. U. Khan, “Joint estimation of amplitude, direction of arrival and range of near field sources using memetic computing,” Progress In Electromagnetics Research C, vol. 31, pp. 199–213, 2012.
- H.-P. Schwefel, Evolution and Optimum Seeking, Sixth-Generation Computer Technology Series, John Wiley & Sons, New York, NY, USA, 1995.
- D. B. Fogel, “The advantages of evolutionary computation,” in Proceedings of the Biocomputing and Emergent Computation (BCEC '97), Skovde, Sweden, September 1997.
- J. Kennedy and R. C. Eberhart, “Particle swarm optimization,” in Proceedings of the IEEE International Conference on Neural Networks, vol. 4, pp. 1942–1948, IEEE, Perth, Australia, December 1995.
- F. Zaman, Estimation of direction of arrival for adaptive beamforming [Ph.D. dissertation], 2013.
- S. H. Ling, K. Y. Chan, F. H. Leung, F. Jiang, and H. Nguyen, “Quality and robustness improvement for real world industrial systems using a fuzzy particle swarm optimization,” Engineering Applications of Artificial Intelligence, 2015.
- R. Storn and K. Price, “Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces,” Journal of Global Optimization, vol. 11, no. 4, pp. 341–359, 1997.
- S. W. Sloan, “A steepest edge active set algorithm for solving sparse linear programming problems,” International Journal for Numerical Methods in Engineering, vol. 26, no. 12, pp. 2671–2685, 1988.
- W. W. Hager and H. Zhang, “A new active set algorithm for box constrained optimization,” SIAM Journal on Optimization, vol. 17, no. 2, pp. 526–557, 2006.
- P. E. Gill and E. Wong, “User's guide for SNCTRL, an interface for SNOPT 7,” Numerical Analysis Report 11-1, Department of Mathematics, University of California, San Diego, Calif, USA, 2011.
- P. E. Gill, W. Murray, and M. A. Saunders, “User's guide for QPOPT 1.0: a fortran package for quadratic programming,” Report SOL 95-4, Department of Operations Research, Stanford University, Palo Alto, Calif, USA, 1995.
- P. E. Gill, W. Murray, and M. A. Saunders, “User's guide for SQOPT version 7: software for large-scale linear and quadratic programming,” Numerical Analysis Report 06-1, Department of Mathematics, University of California, San Diego, Calif, USA, 2006.
- N. I. M. Gould, D. Orban, and P. L. Toint, “GALAHAD, a library of thread-safe fortran 90 packages for large-scale nonlinear optimization,” ACM Transactions on Mathematical Software, vol. 29, no. 4, pp. 353–372, 2003.
- M. A. Z. Raja, “Stochastic numerical treatment for solving Troesch's problem,” Information Sciences, vol. 279, pp. 860–873, 2014.
Copyright © 2015 Fawad Zaman et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.