Codesign of Beam Pattern and Sparse Frequency Waveforms for MIMO Radar
Multiple-input multiple-output (MIMO) radar takes the advantages of high degrees of freedom for beam pattern design and waveform optimization, because each antenna in centralized MIMO radar system can transmit different signal waveforms. When continuous band is divided into several pieces, sparse frequency radar waveforms play an important role due to the special pattern of the sparse spectrum. In this paper, we start from the covariance matrix of the transmitted waveform and extend the concept of sparse frequency design to the study of MIMO radar beam pattern. With this idea in mind, we first solve the problem of semidefinite constraint by optimization tools and get the desired covariance matrix of the ideal beam pattern. Then, we use the acquired covariance matrix and generalize the objective function by adding the constraint of both constant modulus of the signals and corresponding spectrum. Finally, we solve the objective function by the cyclic algorithm and obtain the sparse frequency MIMO radar waveforms with desired beam pattern. The simulation results verify the effectiveness of this method.
Compared with traditional single-station radar and phased array radar, multiple-input multiple-output (MIMO) radar has many advantages: higher resolution to targets; easier detection of the properties of moving targets; better ability to identify parameters; higher degrees of freedom for waveform design. MIMO radar can be classified into two types according to the sparseness of antennas, distributed MIMO radar and centralized MIMO radar, where the main difference between them is as follows: the distributed MIMO radar has larger space between two transceiver antennas; it can provide space diversity gain to improve the detection of scintillating target; the centralized MIMO radar has smaller space between the two transceiver antennas, where each antenna can transmit a different signal waveform, so it has a good waveform diversity capability. The main work of this paper lies in the joint optimization problem of the beam pattern and the radar waveform design for centralized MIMO radar.
In the field of radar, beam pattern design is one of the hot topics in recent years [1, 2]. For the problem of beam pattern design of centralized MIMO radar, linear array model was built in , where a method of designing beam pattern was proposed according to covariance matrix of the transmitted signal by MIMO radar . The method was further improved in  with gradient search algorithm. However, the objective function of this method is in the integral form such that the complexity of solving the problem is very high. In [5, 6] they proposed and obtained the semidefinite quadratic programming expression of the objective function through transforming the form of covariance matrix of the transmitted waveforms. Then, MIMO radar beam pattern synthesis was implemented by semidefinite programming, where the SEDUMI  optimization tool was adopted. In  they studied when desired covariance matrix is known how to design radar signal waveforms such that the synthetic covariance matrix is as close as possible to the desired covariance matrix . The idea in  was to change the constant modulus constraint of radar signals into the peak-to-average-power ratio (PAR) constraint [9, 10]. In this way, constant modulus radar waveforms can be obtained easily. Similarly, a direct method was proposed in  to get the constant modulus signals, where the quasi-Newton method was used to solve the optimization problem. Furthermore, in the in-depth investigations of beam pattern synthesis, more considerations should be counted, such as the power distribution in space following some certain constraints (e.g., beam pattern sidelobe suppression, etc.), and also the synthesis waveforms might meet some practical constraint limits (e.g., sidelobe of autocorrelation function and cross correlation function, etc.). Moreover, the beam pattern design of wideband MIMO radar was discussed in .
The above researches had sufficient studies about the spatial distribution of MIMO radar beam pattern, where the constraints of the main lobe error, the sidelobe error, and the main lobe fluctuation were taken into account. But the spectrum of radar waveforms was not considered. In practical applications of radar and communication, the design of the spectrum of radar waveforms is important. For example, if the spectrum occupation was not considered, it might cause serious interference to civil radio frequency band, since the radar signal is not allowed in some bands (such as the bands for navigation). In this case, radar waveforms which occupy the whole continuous band are not permitted, so the sparse frequency waveforms with several stopbands should be incorporated necessarily . Sparse frequency waveform is a kind of waveforms in which there are several separate stopbands embedded in a continuous band. In the crowded occupied frequency bands such as HF, VHF, and UHF, the application of sparse frequency waveforms has a significant practical importance. The design method of sparse frequency radar waveforms was studied in [13–16], where power spectra density matching method was proposed to design sparse frequency waveforms with constraint on the sidelobe of the autocorrelation function.
In this paper, we extend the application of sparse frequency waveform to beam pattern design of MIMO radar. First, we calculate the covariance matrix with semidefinite constraint as in . Then, we generalize the objective function of sparse frequency radar waveforms with the known covariance matrix and solve the optimization problem by changing the constant modulus constraint into PAR constraint as in . After that, we design the waveforms with sparse frequency spectrum by power spectra density matching method with the quasi-Newton method. At last, we use cyclic algorithm to obtain sparse frequency radar waveforms with desired beam pattern. In simulations, we compare the acquired beam pattern with the proposed method to the ideal beam pattern and the beam pattern with desired covariance matrix , which show that they are very close to each other. Then we simulate the power spectra density and autocorrelation function of the acquired sparse frequency radar waveforms and give the corresponding analysis of the properties of spectrum and autocorrelation function.
2. MIMO Radar Signal Model
Consider a centralized MIMO radar system with radar transmitters. The number of samples in a transmitted signal pulse repetition period is . Let be a discrete-time radar transmitted baseband signal. The transmitted signal of the th radar transmitter isand a transmitted signal pulse of the MIMO radar system iswhere the superscript symbol denotes the transpose. The matrix with size represents MIMO radar waveforms matrix; that is,According to (1)–(3), the relationship among , , and is
The wavelength of the signal transmitted in the MIMO radar system is , where the distance between each of the two radar transmitters is . The linear array of radar transmitters is shown in Figure 1.
As the uniform linear array (ULA) shown in Figure 1, the distance between the two transmitters is a half wavelength. Let represent the angle parameter of the target of interest. Then the signal at the target can be expressed aswhere represents radar carrier frequency and denotes the time delay of the th antenna transmitted signal to the target. If there is no change of the relative location between the known target and MIMO radar array, the time delay is a constant. We introduce a symbol to represent both the radar carrier frequency and the time delay, which isWhen the radar carrier frequency is known, is dependent only on . Then from (5) and (6), the interested signal can be expressed aswhere the superscript symbol denotes the conjugate transpose. It can be seen from (7) that, without the attenuation caused by distance, the power radiated at the target by radar array in a pulse repetition period should bewhere is the covariance matrix of the transmitted signal , which is
We can observe that is a function of from (8), which means it actually reflects the radar transmitting energy distribution in space. Hence, it is called “transmitting beampattern.” According to (8) and (9), the problem of how to design MIMO radar transmitting beam pattern can be converted to the problem of how to design the covariance matrix .
The value of covariance matrix has a direct connection to the performance of the radar system. For example, when the rank of is “1,” the transmitted signals are perfectly correlated, which is equivalent to the case of conventional phased array radar; when is a diagonal matrix, the transmitted signals are mutually orthogonal. In a practical MIMO radar system, the transmitted signals are usually not mutually orthogonal, but there is a certain correlation between them. Through changing matrix , we can design the beam pattern with different properties, depending on the requirements.
Given the ideal power distribution , we can design the covariance matrix such that the transmitting beam pattern is close to the ideal beam pattern. The optimization of this problem  iswhere the optimization variable is a ratio which relates the magnitude of the beampattern. means the equal power constraint of each transmitter channel, and means the semidefinite constraint of covariance matrix. It can be seen that (10) is a problem of semidefinite quadratic programming. We can solve it by SEDUMI tool to get the covariance matrix such that it can be close to the ideal beam pattern.
Since the covariance matrix is assumed to be known, the beam pattern of the radar waveforms with constraints would not be closer to the ideal beam pattern than . Thus, we call the “optimal covariance matrix.”
The main topic of this paper is the codesign of beam pattern and sparse frequency waveforms for MIMO radar. We start from the optimal covariance matrix and then make the acquired beam pattern close to the desired beam pattern and get the radar waveforms with certain constraints (constant modulus constraint and sparse frequency constraint). The problem of how to get the waveforms with constant modulus constraint was discussed in . We first introduce briefly the method of obtaining constant modulus signal and then propose the method of imposing sparse frequency constraint, that is, power spectra density matching.
3. Sparse Frequency Radar Waveforms Design
3.1. Constant Modulus Waveforms Design
First we analyze the relationship between the covariance matrix and the radar waveforms matrix . From (9), we haveWhen the covariance matrix is known and there is no constraint, and are one-to-many relationship. In general, the sample number of radar waveform pulse is greater than the number of radar transmitters; that is, , so is a tall matrix. The problem is how to calculate the tall matrix from the square matrix . We introduce an matrix , which satisfies . Also, the covariance matrix has the same dimension as the unit matrix. From and (11), we have
Based on (12) we can observe that when the covariance matrix is known we can find a way to calculate the radar waveforms matrix through the converting matrix . At the same time it is worth noting that the choice of matrix will decide the properties of radar waveforms matrix . In a practical radar system, the main concern is the constant modulus of radar waveform signals. Let represent the set of that meet the constraint of constant modulus in practice. Then the optimization problem that calculates the radar waveforms matrix through covariance matrix should bewhere is the unit matrix constraint of matrix and is the constant modulus constraint of matrix . Objective function (13) can be solved by cyclic algorithm. The basic idea is as follows.(1)Fix the matrix and minimize in terms of the matrix .(2)Fix the matrix and minimize in terms of the matrix .
These two steps should be implemented alternately.
Stop criteria are as follows: check the Frobenius norm value which is the difference of matrix in two adjacent iterations. When is smaller than the threshold, stop the iteration. The radar waveforms matrix with constant modulus can be obtained. The specific steps will be discussed in Section 3.2.
3.2. Sparse Frequency Waveforms Design
How to get a waveform sequence with sparse frequency spectrum from a constant modulus sequence? We can use a method of power spectra density matching to realize it. For the continuous signal process, we assume that represents a radar waveform, where the duration is . Given an ideal power spectra density with sparse frequency spectrum, expressed by , the optimization problem of designing sparse frequency waveform through power spectra density matching can be expressed aswhere is the Fourier transform operation of and is the constant modulus constraint of signal. Note that it is the continuous form of signal waveform.
In order to get the corresponding discrete representation, we use sequence to represent a sparse frequency radar waveform. Each element in could be written aswhere the effect of is to normalize waveform energy and is the phase value of radar waveform at the th moment. We can use to represent the phase vector composed by phase values; that is,Since every element in is only related to the phase in the current moment, the radar waveform can be written as a function of . Let be the discrete representation form of ideal power spectra density . Then the objective function of solving sparse frequency waveform with known can be expressed bywhere the superscript symbol is the conjugate operation, means the Hadamard multiplication between two matrices, and is the discrete Fourier transform matrix. Every element in is defined as , .
Since the objective function in (17) is a quartic unconstrained nonconvex optimization problem, it will be difficult to get global optimal solution with general approach. Thus, we can get the local optimal solution by adopting the quasi-Newton method. Compared with the Newton method, quasi-Newton method has superlinear convergence rate with less computational complexity. Therefore, besides quadratic convex optimization problem, it is widely used for general optimization problem. Similar to steepest descent method, the quasi-Newton method is also required to calculate the gradient of the objective function in each iteration. In the following discussion, the objective function is transformed first, and then quasi-Newton method and BFGS correction will be used to solve the problem.
Let represent the objective function and , . Equation (17) can be written asSince is known, we can combine and . Letting , then we haveThen the gradient of the objective function on can be calculated bywhere the symbol means the operation of constructing a diagonal matrix and the symbol means the imaginary part of a complex number.
We need the gradient value of the objective function at each iteration. Let subscript symbol represent the index of each iteration. Then the steps for solving the objective function are as follows.
Step A. Initialize and a tolerance value . Set and . Calculate by (20).
Step B. Get the search direction , calculate , and use linear search method to find the step size which can minimize . Let and .
Step C. If , then stop the iteration, and output the result ; if , then go to Step D.
Step D. Use (20) to calculate ; let ; then calculate , which is
Step E. Let ; go back to Step B.
During the process, is used for comparison to the Frobenius norm of , that is, the difference of matrix in two adjacent iterations, which is used as the condition of the iteration termination; is the approximation of the inverse of the Hessian matrix; the calculation of is consistent with BFGS correction, and more issues in detail are referred to in [17, 18]. From the above steps, the main computation load of quasi-Newton method and BFGS correction lies in the computation of the objective function as well as the . After iterative calculation of waveform with initial phase vector , we can get the radar waveform with result phase vector . has the sparse frequency spectrum. Because the iteration process only changes the value of the phase vector , there is no effect on the modulus value of the signal. Therefore, if the initialization sequence is constant modulus, then the acquired waveform still keeps the constant modulus.
The number of transmitters of MIMO radar system is , so we should require that each transmitter has the same constraint of sparse frequency. Set the phase vector of waveform of the th radar transmitter to be . We get the transmitted waveform as . We assume that radar transmitters have the same constraint of sparse frequency, meaning that the corresponding ideal power spectra density is . Then the objective function of sparse frequency waveforms design of MIMO radar system becomeswhere mean that the optimization variables of the objective function (22) are the phase vectors of radar transmitted waveforms.
After adding sparse frequency condition, let the th transmitted waveform of MIMO radar system transmitters be . The objective function of constant modulus waveforms with the sparse frequency spectrum obtained from the optimal covariance matrix can be written as
We deal with the objective function (23) with cyclic algorithm by introducing a value as the threshold of the iteration termination. Then the solving steps can be written as follows.
Step 1. Initialize matrix and tolerance value .
Step 2. Fix as known; calculate the matrix to minimize , where the constraint is satisfied.
Step 3. Fix as known; calculate the matrix to minimize , where the constant modulus constraint is satisfied.
Step 4. Let , and calculate the phase vectors of radar waveforms. With the help of solving sparse frequency waveform in (20), we can get the optimal phase vectors and update phase vectors and update according to the formula .
Step 5. Check the Frobenius norm value of the difference of matrix in two adjacent iterations. If , stop the iteration and output the result; if , then go back to Step 2.
In Step 2, the purpose is to obtain the needed matrix from matrix . Let the estimated value of be . If the matrix is given, we havewhere is a constant, symbol denotes the trace of a matrix, and symbol means the real part of a complex number. Consider in (24) as a whole, and letwhich means the singular value decomposition (SVD) of . We can have the special structure of this SVD decomposition, where is a diagonal matrix, is the matrix after SVD decomposition, and is the matrix after SVD decomposition (e.g., see ). Then in (24) can be expressed in terms of , , and ; that is,where denotes the diagonal element of matrix and means the diagonal element of . Since has the inequality relationFrom (24)–(27), we can get the results as follows:For , we have the inequalityBased on the upper bound of (29), the estimated value of isDuring Step 2, we can get the needed matrix that holds the orthogonal constraint.
In Step 3, the goal is to obtain the needed matrix from matrix . Let the estimated value of be . Since the matrix obtained from Step 2 is known, the matrix isFor example,does not meet the requirement of constant modulus. We could adjust its amplitude by the following definition:where is a constraint with range , symbol means the angle of a vector, and is the weighted value which is expressed bywhere is the set that represents least magnitude entry in , the number of elements in is , is a variable which belongs to , and the value of is decided by (e.g., see  for details).
From (33) and (34), we obtain the new waveform through adjusting the amplitude of waveform . Here, the idea is to convert the constant modulus constraint to the peak-to-average-power ratio (PAR) constraint, which isIf is a small value, the acquired waveform after the optimization is close to constant modulus. If , the constraint is equivalent to the constant modulus constraint, in some modern radar systems (e.g., see ). radar waveforms in MIMO radar system can be calculated, respectively, by the PAR constraint. Hence, the constant modulus constraint of MIMO radar waveforms can be written as PAR constraints:
In Step 4, the purpose is to obtain a new waveform with sparse frequency spectrum from the acquired constant modulus waveforms matrix . Let the estimated value of sparse frequency waveform be . Consider the constant modulus waveforms matrix obtained from Step 3 as known. Then we have the relationshipto obtain phase vectors . We set the value from formula (37) as the initial value and start from . The algorithm is similar to the single sparse frequency waveform. The steps are as follows.
Step a. Initialize and tolerance value , set , and calculate with (20).
Step b. Let the search direction , calculate , and use linear search method to find the step size which can minimize ; let and let .
Step c. If , output the result and go to Step f; if , go to Step d.
Step d. Calculate with formula (20), let , and calculate :
Step e. Let and go back to Step b.
Step f. If , then stop the iteration; if , let ; go back to Step a.
Compared with the process for the single sparse frequency waveform, solving MIMO radar waveforms takes the original iterative calculation as a submodule, which actually just adds the outer loop by increasing variable . This process is to get the optimal phase vectors and the corresponding waveforms . Obviously, the steps written above are described in a sequential processing order. If multiple processors are allowed, then the optimization of radar waveforms can be performed in parallel. It is worth noting that in Step 4 we only change the value of waveforms matrix through phase vector matrixes , so the constant modulus property in Step 3 still holds.
After Step 4, we obtain sparse frequency waveforms. Then, we can calculate the Frobenius norm value of the difference of matrix in two adjacent iterations and use it as a termination condition.
If , output as the optimal result. Letting the covariance matrix of be , we can get the sparse frequency radar waveforms with constant modulus.
In this way, the beam pattern of its covariance matrix is very close to the beam pattern of the optimal covariance matrix .
3.3. Chart of the Design Steps
Here we use the chart to show the basic steps of codesign of the transmitted beam pattern and sparse frequency radar waveforms proposed in this paper.
As shown in Figure 2, first we use the semidefinite quadratic programming tool SEDUMI to solve the optimal covariance matrix . Then we start from optimal covariance matrix and introduce matrix and ideal power spectra density to minimize the objective function. At last we solve the optimization problem by cyclic algorithm, which solves the unit matrix constraint of matrix , constant modulus constraint of matrix , and sparse frequency constraint. Sparse frequency constraint of waveform , , can be solved by power spectra density matching with the quasi-Newton method, as shown in the dashed box of Figure 2. In the step of sparse frequency constraint, the optimization of radar transmitted waveforms can be calculated, respectively. The iteration stops if termination condition is satisfied, and then output the waveform matrix .
In this section, we numerically verify the effectiveness for the proposed codesign method of MIMO radar beam pattern matching and sparse frequency radar waveforms. Above all, we find the optimal covariance matrix by optimization tool according to the literature . Second, we calculate the sparse frequency waveforms from the optimal covariance matrix with the method proposed in this paper. Finally, we compare the beam pattern of covariance matrix of the calculation results to the one with optimal covariance matrix .
Assume that the MIMO system has radar transmitters and the length of each transmitter transmit waveform sequence is . We set the desired ideal beam pattern through angle function. The angle function iswhere is the central position of the desired ideal beam pattern, here it is assumed to be , meaning that the central position is , and is the half beam-width of the desired ideal beam pattern, here it is , meaning that the beam-width is . Substitute into formula (10) and solve the problem with the SEDUMI tool . The optimal covariance matrix can be obtained, and the corresponding beam pattern is shown in Figure 3.
Figure 3 shows the beam pattern of the variance matrix by solving semidefinite quadratic programming problem with SEDUMI optimization tool. The horizontal axis represents the angle and the vertical axis represents the normalized logarithm of the beam pattern power. From the figure we can observe that the desired ideal beam pattern is wide beam pattern with central position and beam-width 60°. The beam pattern of the optimal covariance matrix by semidefinite quadratic programming has a strong power distribution in space −30°~+30° section but little power in other sections. The power of its peak sidelobe is lower than −10 dB. The beam pattern of has smooth main lobe, which shows good agreement with the desired ideal beam pattern.
The steps of solving sparse frequency waveforms from covariance matrix are as follows. Set a tolerance value and ; the power spectra density of ideal sparse frequency is assumed to have 3 passbands and 2 stopbands. The gain in passband is 5 dB and in stopband is −20 dB; that is, where symbol is the set of passbands. Because the number of waveform samples is , we can set . We have an ideal power spectra density of the narrowband radar waveform which is shown in Figure 4.
Figure 4 illustrates the desired ideal power spectra density set. We can observe that the spectrum of the radar waveform has 3 passbands and 2 stopbands, which is the goal of designing sparse frequency waveforms. The initial of can be realized by random function. Through implementing the steps by (23) in Section 3.2, we can get the result which represents the waveforms from radar transmitters.
Figure 5 shows the power spectra density of the obtained radar waveforms in frequency domain. In Figure 5(a), coordinate represents the index of the transmitters in MIMO radar system, where in this experiment there are 10 transmitters; coordinate represents the frequency with unit MHz; coordinate shows the amplitude of the power spectra density with logarithmic representation. In order to facilitate analysis, Figure 5(b) combines these 10 curves of radar waveforms in projection form, so we could observe the common properties. It can be seen from the figure that the waveforms calculated by cyclic algorithm have the sparse frequency property. The spectrum has 3 passbands and 2 stopbands with good shape at the boundaries. Also, the 10 power spectra density curves have similar shapes. From the practical simulation data, we can calculate that the average value of the 10 power spectra density in the passbands is 5.84 dB, and the average value in the stopbands is −20.63 dB, which can meet the requirement. After obtaining sparse frequency waveforms, we obtain the covariance matrix from waveforms matrix .
Figure 6 shows the comparison of beam pattern of resulting sparse frequency waveforms, ideal beam pattern, and beam pattern of optimal covariance matrix , where the horizontal axis is the angle and the vertical axis is the normalized logarithm of the beam pattern power. It can be observed from Figure 6 that the resulting beam pattern is very close to the beam pattern with optimal covariance matrix in the region of high power, where the peak sidelobe of the beam pattern with optimal covariance matrix is −10.2 dB and the peak sidelobe of the beam pattern with sparse frequency waveforms by the proposed method is −9.3 dB. It can be concluded that MIMO radar transmitted waveforms with sparse frequency properties can be obtained successfully by the method proposed in this paper and the corresponding covariance matrix is close to the optimal covariance matrix .
Next, we analyze the autocorrelation function of the radar waveforms .
Figure 7 shows autocorrelation function curves of sparse frequency waveforms . In Figure 7(a), coordinate represents the index of the transmitters in MIMO radar system; coordinate denotes the time delay with normalized unit; coordinate means the magnitude of the autocorrelation function with logarithmic representation. In Figure 7(b) we project these autocorrelation function curves on the plane such that we can observe the common properties. It can be seen from the figure that these autocorrelation function curves have the same trend. According to the practical simulation data, we calculate the average peak sidelobe of this autocorrelation function, which is −20.3 dB. It can be concluded that the resulting sparse frequency radar waveforms show good autocorrelation properties.
The simulation above is only for single wide beam pattern. For the case of multibeam, we can have similar results. We assume that there are 3 beams of the desired ideal beam pattern, with central positions , , and , respectively, and with the same half beam-width . We can perform the simulation and obtain results as shown in Figure 8.
Figure 8 shows, in the case of multibeam, the comparison of beam pattern with sparse frequency waveforms, ideal beam pattern, and the beam pattern with optimal covariance matrix . It can be apparently seen that the resulting beam pattern is still close to the beam pattern of optimal covariance matrix calculated by the SEDUMI tool. Therefore, the effectiveness of the method proposed for the multibeam case can be verified.
This paper studied the codesign of beam pattern and sparse frequency radar waveforms for centralized MIMO radar system. We started from the optimal covariance matrix , introduced the constraint of sparse frequency, and considered the constant modulus constraint of practical radar waveforms. Then we minimized the objective function by cyclic algorithm. At last we got the sparse frequency radar waveforms with constant modulus. The beam pattern of covariance matrix of the radar waveforms is close to the optimal covariance matrix . The spectrum of the radar waveforms with the proposed method has the sparse property, which can satisfy the specified requirement. Simulation results verify that the proposed method can provide good properties of beam pattern with good sparse spectrum property, where the acquired radar waveforms also show good autocorrelation property.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
This work was supported by the National Natural Science Foundation of China (Grant no. 61471019).
D. R. Fuhrmann and G. san Antonio, “Transmit beamforming for MIMO radar systems using partial signal correlation,” in Proceedings of the 38th Asilomar Conference on Signals, Systems and Computers, pp. 295–299, Pacific Grove, Calif, USA, November 2004.View at: Google Scholar
J. A. Tropp, I. S. Dhillon, R. W. Heath Jr., and T. Strohmer, “CDMA signature sequences with low peak-to-average-power ratio via alternating projection,” in Proceedings of the 37th Asilomar Conference on Signals, Systems and Computers, vol. 1, pp. 475–479, Pacific Grove, Calif, USA, November 2003.View at: Google Scholar
G. H. Wang and Y. L. Lu, “Designing sparse frequency waveform using iterative algorithm,” in Proceedings of the 11th International Radar Symposium (IRS '10), pp. 1–4, IEEE, Vilnius, Lithuania, June 2010.View at: Google Scholar
J. Nocedal and S. J. Wright, Numerical Optimization, Springer, 2nd edition, 2006.
A. Antoniou and W. S. Lu, Practical Optimization: Algorithms and Engineering Applications, Springer, 2007.