序列二次规划算法

1. Legendre PolynomialsNamed in honor of Adrien-Marie (1752-1833) the mathematician, not Louis (1752-1797) the politician.1.1. IntroductionA polynomial is a finite sum of terms like akxk, where k is a positive integer or zero. There are sets of polynomials such that the product of any two different ones, multiplied by a function w(x) called a weight function and integrated over a certain interval, vanishes. Such a set is called a set of orthogonal polynomials. Among other things, this property makes it possible to expand an arbitrary function f(x) as a sum of the polynomials, each multiplied by a coefficient c(k), which is easily and uniquely determined by integration. A Fourier series is similar, but the orthogonal functions are not polynomials. These functions can also be used to specify basis states in quantum mechanics, which must be orthogonal.1.2. Legendre Polynomials DenifitionThe Legendre polynomials Pn(x), n = 0, 1, 2 . are orthogonal on the interval from -1 to +1, which is expressed by the integral.The Kronecker delta is zero if n m, and unity if n = m. In most applications, x = cos , and varies from 0 to . In this case, dx = sin d, of course. The Legendre polynomials are a special case of the more general Jacobi polynomials P(,)n(x) orthogonal on (-1,1). By a suitable change of variable, the range can be changed from (-1,1) to an arbitrary (a,b). The weight function w(x) of the Legendre polynomials is unity, and this is what distinguishes them from the others and determines them.1.3. ApplicationsThe Lengendre polynomials are very clearly motivated by a problem that often appears. For example, suppose we have an electric charge q at point Q in the figure at the left, one of a group whose positions are referred to an origin at O, and we desire the potential at some point P. The distance PO is taken as unity for convenience; simply multiply all distances by the actual distance PO in any particular case. The potential due to this charge is q/R. We can find R as a function of r and by the Law of Cosines: R2 = 1 + r2 - 2r cos = 1 - 2rx + r2, where x = cos . Now we expand 1/R in powers of r, finding 1/R = Pn(x)rn. The function 1/R is called the generating function of the Legendre polynomials, and can be used to investigate their properties. Generating functions are available for most orthogonal polynomials, but only in the Legendre case does the generating function have a clear and simple meaning.If we let x = 1, we find that Pn(1) = 1, and Pn(-1) = (-1)n. By taking partial derivatives of 1/R with respect to x and r, and then considering the coefficients of individual powers of r, we can find a number of relations between the polynomials and their derivatives. These can be manipulated to find the recursion relation, (n + 1) Pn+1(x) = (2n + 1)x Pn(x) - n Pn-1(x), and the differential equation satisfied by the polynomials, (1 - x2) Pn(x) -2x Pn(x) + n(n + 1) Pn(x) = 0. The recurrence relation allows us to find all the polynomials, since it is easy to find that P0(x) = 1, P1(x) = x directly from the generating function, and this starts us off. The differential equation allows us to apply the polynomials to problems arising in mathematics and physics, among which is the important problem of the solution of Laplaces equation and spherical harmonics.The recurrence relation shows that the coefficient An of the highest power of x satisfies the relation An+1 = (2k + 1)/(k + 1) An, and so from the known coefficients for n = 0, 1 we can find that the coefficient of the highest power of x in Pn is 1.3.5.(2n-1)/n!.The polynomials can also be found by solving the differential equation by determining the coefficients of a power series substituted in the equation. This method was often used in quantum mechanics texts (see Reference 3), since the students were not usually acquainted with the mathematics of orthogonal polynomials. This method does not allow one to investigate the properties of the polynomials in any detail, however, yielding only the individual polynomials themselves.Consider the polynomials Gn(x) = dn/dxn (x2 - 1)n. The quantity to be differentiated is indeed a polynomial, of degree 2n, and consisting of only even powers. When differentiated n times, it becomes a polynomial of order n consisting of either all odd or all even powers of x, as n is odd or even. The coefficient of the highest power of x is 2n(2n-1)(2n-2).(n+1), and the first two polynomials are 1 and 2x. If G(x) is substituted in the recurrence relation for the Legendre polynomials, it is found to satisfy it. If we divide G(x) by the constant 2nn!, then the first two polynomials are 1 and x. Therefore, Pn(x) = (1/2nn!) dn/dxn (x2 - 1)n. This is called Rodriguess formula; similar formulas exist for other orthogonal polynomials.The great advantage of Rodrigues formula is its form as an nth derivative. This means that in an integral, it can be used repeatedly in an integration by parts to evaluate the integral. The orthogonality of the Legendre polynomials follows very quickly when Rodrigues formula is used. There is a Rodrigues formula for many, but not all, orthogonal polynomials. It can be used to find the recurrence relation, the differential equation, and many other properties.For finding solutions to Laplaces equation in spherical coordinates, the Legendre polynomials are sufficient so long as the problem is axially symmetric, in which there is no -dependence. The more general problem requires the introduction of related functions called the associated Legendre functions that are actually built up from Jacobi polynomials, and can also be expressed in terms of derivatives of the Legendre polynomials. Physics texts generally approached the problem from first principles, never mentioning Jacobi polynomials, and thereby losing valuable insight.The Jacobi polynomials P(,)n(x) are orthogonal on (-1,1) with weight function w(x) = (1 - x)(1 + x). Their Rodrigues formula is P(,)n(x) = (-1)n/2nn! (1 - x)-(1 + x)- dn/dxn (1 - x)+n(1 + x)+n. The ordinary Legendre polynomial Pn(x) = P(0,0)n(x). They satisfy the differential equation (1 - x2)P(,)n + - - ( + + 2)x P(,)n + n( + + n + 1) P(,)n = 0.In solving Laplaces equation by the method of separation of variables, one obtains for the dependence T(x), x = cos , the differential equationd/dx(1 - x2)dT/dx = l(l+1) - m2/(1 - x2T = 0The substitution T(x) = (1 - x2)m/2y(x) now gives the equation(1-x2)y - 2(m + 1)xy + l(l+1) - m(m+1)y = 0,which we recognize as satisfied by the Jacobi polynomial P(m,m)l-m(x). Hence, T(x) = (1 - x2)m/2y(x) P(m,m)l-m(x). This is the associated Legendre function, often denoted Pml(x) in physics texts (e.g., Reference 4), and defined there as (-1)m(1 - x2)m/2 dm/dxm Pl. The subscript is no longer the degree of the polynomial.All the above is for a positive m. Since the equation contains m2, the solution for negative m is essentially the same, except perhaps for a multiplicative factor. This is of little consequence for the traditional applications of spherical harmonics, but is critical for quantum mechanics, where relative phases matter. The choice in physics is that P-ml(x) = (-1)m(l - m)!/(l + m)! Pml(x), where m is always positive on the right. If you work the functions out explicitly, you will find that the functions for +m and -m are essentially the same, as might be expected, and differ at most by a factor of -1.For the same m, Pml(x) and Pml(x) are orthogonal, and the integral of the square of Pml(x) is the same as for Pl(x), multiplied by (l - m)!/(l + m)!. The functions are not orthogonal for different values of m; orthogonality of the spherical harmonics in this case depends on the functions.1.4. References1. M. Abramowitz and I. Stegun, Handbook of Mathematical Functions (Washington, D.C.: National Bureau of Standards, Applied Mathematics Series 55, June 1964). Chapter 22.2. D. Jackson, Fourier Series and Orthogonal Functions (Mathematical Assoc. of America, Carus Mathematical Monographs No. 6, 1941). Chapter X.3. L. Pauling and E. B. Wilson, Introduction to Quantum Mechanics (New York: McGraw-Hill, 1935). Chapter V.4. J. D. Jackson, Classical Electrodynamics, 2nd . ed. (New York: McGraw-Hill, 1975), Chapter III.2. Gauss型积分2.1. Gauss型求积公式的构造方法(1)求出区间a,b上权函数为W(x)的正交多项式pn(x) (2)求出pn(x)的n个零点x1 , x2 , xn 即为Gsuss点. (3)计算积分系数 2.2. 几种Gauss型求积公式2.2.1. Gauss-Legendre求积公式区间-1,1上权函数W(x)=1的Gauss型求积公式，称为Gauss-Legendre求积公式，其Gauss点为Legendre多项式的零点。公式的Gauss点和求积系数可在数学用表中查到。由于因此，a,b上权函数W(x)=1的Gauss型求积公式为2.2.2. Gauss 公式的余项2.3. Gauss-Laguerre求积公式区间0,+)上权函数W(x)=e-x的Gauss型求积公式，称为Gauss-Laguerre求积公式，其Gauss点为Laguerre多项式的零点。公式的Gauss点和求积系数可在数学用表中查到。由，所以，对0, +)上权函数W(x)=1的积分，也可以构造类似的Gauss-Laguerre求积公式：。2.4. Gauss-Hermite求积公式区间（-, +)上权函数W(x)=e-x2的Gauss型求积公式，称为Gauss-Hermite求积公式，其Gauss点为Laguerre多项式的零点。公式的Gauss点和求积系数可在数学用表中查到。3. Legendre多项式3.1.1. 隐式表达式3.1.2. 显式表达式其中正交性质：递推关系：3.2. 高斯积分的计算过程（1）、将a,b区间积分变换到-1,1区间；（2）、求高斯点及权wi(i=1,2,n)；（3）、代入高斯公式：，其中，最终，4. 高斯伪谱法的轨迹优化方法下面讨论基于高斯伪谱法的质点弹道优化方法。4.1. 质点运动方程取，则考虑地球自转，忽略地球扁率，不考虑侧滑情况下的三自由度质点运动模型可以描述为4.2. 约束条件，具体介绍如下4.2.1. 过程约束热流、过载、动压以及拟平衡滑翔条件等。其中动压、热流和过载必须严格满足。4.2.1.1. 热流密度研究飞行器轨迹优化时通常以驻点热流作为约束条件，因为驻点是飞行器加热较为严重的区域。其工程估算表达式为：对于高超声速飞行器，一般取，。是飞行器头部半径相关的参数。按照常用的热流密度表达式，热流密度约束为：为飞行器头部曲率半径，与飞行器特性相关的参数。4.2.2. 动压约束考虑到动压对飞行器控制系统的影响和侧向稳定性的要求，再入过程中动压满足4.2.3. 过载约束为了保证结构安全，需考虑过载约束。弹箭类飞行器一般对法向过载进行约束。而升力体飞行器其法向和轴向都可能产生大过载。因此对总过载进行约束：4.2.3.1. 拟平稳滑翔条件相对于前面几种约束，拟平稳滑翔约束不太严格，是一种考虑飞行器控制能力的“软约束”，即所谓的“no-skip”条件，保证飞行器不再跳跃出大气层。无因次化后的表达式为：其中是考虑哥氏加速度和牵连加速度的附加项，一般再入问题中其取值较小；是平衡滑翔边界对应的倾侧角。考虑，则上述约束简化为4.2.3.2. 控制变量约束4.2.4. 终端约束条件4.2.4.1. 到达指定点的约束，对于对敌打击武器平台，还附加终端航迹角和速度的约束，4.2.4.2. 到达指定区域的约束定义飞行器距离目标点的地面航程为，则约束条件为，定义航向角与相对目标点的视线角的偏差为，为使末制导段能更准确的到达目标，要求到达指定区域时航向也满足一定约束，即4.2.5. 目标函数可以是下面的一个或多个的加权和。4.2.5.1. 吸热最少4.2.5.2. 弹道平滑性4.2.5.3. 射程最大4.2.6. 轨迹优化数学模型状态变量为，。4.2.7. 轨迹优化问题的转化4.2.7.1. 区间变换将原轨迹优化问题的时间区间变换到GPM求解需要的区间。因此令，优化问题化为4.2.7.2. 状态变量和控制变量的多项式近似设利用高斯伪谱法对优化问题进行离散处理选取的插值点个数为，分别为，其中，。（）为阶Legendre正交多项式的零点，这些点位于。通过上述定义，则状态变量和控制变量的拉格朗日多项式近似为其中和分别为以和为节点的拉格朗日插值基函数：而终端状态约束转化为：其中为高斯加权系数：在状态变量、控制变量用拉格朗日多项式近似，终端约束用Gauss-Legendre积分近似后，下一步建立NLP问题。4.2.7.3. NLP问题的建立状态微分用拉格朗日多项式插值结果求导得到：于是用上式代替质点运动学方程的导数，有在插值节点（）处离散化得到：，其中所以于是原问题的微分方程约束转化为多项式约束。将边值约束化为将过程约束转化为，经过上述转化，轨迹优化问题的求解就转化为NLP问题的求解。即求插值点处的状态变量、，控制变量以及初始时刻和终端时刻，使得满足：（1）终端约束：；（2）动力学方程约束：，（3）边值约束：（4）过程约束：，其中（5）性能指标最优：4.2.8. 标准NLP问题一般，因此只有终端时间需要寻优。于是令则轨迹优化问题转换为标准的NLP问题。其标准优化模型为：可见，决策变量共有个；等式约束共，为边值约束方程的个数；不等式约束个数等于过程约束条件的个数。上述问题简化描述为此处，4.2.9. NLP问题的求解目前对NLP问题求解的方法中，SQP算法一起整体的收敛性和局部超一次收敛，被认为是目前求解NLP问题的最好方法之一。此处采用SQP算法对上述问题进行求解。前述NLP问题的拉格朗日函数为利用一阶泰勒级数近似约束条件，二阶泰勒级数近似目标函数，则NLP问题等价于求解下面的QP问题：其中，为拉格朗日函数的Hessian矩阵的近似正定矩阵，采用改进的BFGS公式对其进行更新：，为不大于的正数，。和分别为步长因子和迭代的搜索方向。算法步骤：Step1：给定初始点，初始矩阵及控制误差；令，其中插值点处的可根据给定状态和控制变量的初值及终值，通过插值计算获得；Step2：求解QP问题，确定搜索方向及相应的拉格朗日乘子；Step3：由精确一维搜索确定补偿因子，从而计算出新的迭代点；Step4：若收敛判别满足，则取NLP问题的最优解为，停止计算；否则，转Step5；Step5：更新Hessian矩阵，令，返回Step2。求的NLP问题的最优解后，可得轨迹优化问题状态变量和控制变量在插值点处的最优值及末端时刻。需要进一步考虑的问题：如何给定？；精确一维搜索？；QP问题的求解？函数导数的基于有限差分计算初值选为正定矩阵即可，如选择单位矩阵参见matlab的函数fmincon和nlconst（C:MATLAB6p5toolboxoptimprivate）qpsub