3.2 非线性规划L I NGO程序设计基础

线性规划的应用范围十分广泛，但仍存在较大局限性，对许多实际问题不能处理。非线性规划比线性规划有着更强的适用性。事实上，客观世界中的问题多是非线性的，给予线性处理大多是近似的，是在做了科学的假设和简化后得到的。在实际问题中，有一些不能进行线性化处理，否则将严重影响模型对实际问题近似的可依赖性。但是，非线性规划问题在计算上通常比较困难，理论上的讨论，也不能像线性规划的问题那样，给出简捷的形式和透彻全面的结论。

例3.3 板材切割问题

某钢管零售商从钢管厂进货，将钢管按照顾客的要求切割后售出。从钢管厂进货时得到的都是长度为19m原料钢管。现有一个客户需要长度为4m的钢管50根、长度为6m的钢管20根和长度为8m的钢管15根。应如何下料最节省？

零售商如果采用的不同切割模式太多，将会导致生产过程的复杂性，从而增加生产和管理成本。所以，该零售商规定采用的不同切割模式不能超过3种。此外，该客户除需要上述中的三种钢管外，还需要长度为5m的钢管10根，应如何下料最节省？

模型设计

问题一

首先，应当确定哪些切割模式是可行的，所谓一个切割模式，是指按照客户需要在原料钢管上安排切割的一种组合。例如，我们可以将长度为19m的钢管切割成3根长度为4m的钢管，余料为7m；或者将长度为19m的钢管切割成长度为4m、6m和8m的钢管各1根，余料为1m。显然，可行的切割模式可以有很多。

其次，应当确定哪些切割模式是合理的。通常，假设一个合理的切割模式的余料不应该大于或等于客户需要的钢管的最小尺寸。例如，将19m长的钢管切割成3根4m的钢管是可行的，但余料为7m，可以进一步将7m的余料切割成4m钢管（余料为3m），或者将7m的余料切割成6m钢管（余料为1m）。在这种合理性假设下，切割模式一共有7种，如表3-6所示。

问题化为在满足客户需要的条件下，按照哪些合理的模式，切割多少根原料钢管最为节省。所谓节省可以有两种标准：一是切割后剩余的总余料最小，二是切割原料钢管的总根数最少。下面将对这两个目标分别讨论。

决策变量：用xi表示按照第i种模式（i=1,2, …,7）切割的原料钢管的根数，显然它们应当是非负数。决策目标：以切割后剩余的总余量最小为目标，由表3-6可得到：

minZ1=3x1+x2+3x3+3x4+x5+x6+3x7

以切割原料钢管的总根数最少为目标，则有：

minZ2=x1+x2+x3+x4+x5+x6+x7

约束条件为满足客户的要求，则有：

【思考题】 请思考两种目标函数下，取得的最优方案是否相同？为什么？如果两种最优方案不同，请思考何种更为合理？

问题二

按照解问题一的思路，可以通过枚举法首先确定哪些切割模式是可行的，但由于需求的钢管规格增加到4种，所以枚举法的工作量较大。

通过枚举法可以确定16种可行且合理的切割模式。决策变量：用xi表示按照第i种模式（i=1,2, …,16）切割的原料钢管的根数，显然它们应当是非负数。决策目标：以切割原料钢管的总根数最少为目标，则有：

分析第二问的难点在于：该零售商规定采用的不同切割模式不能超过3种。为了实现该要求，引入额外的决策变量：用ri 表示采用第i种模式（i=1,2, …,16）。

因此，对切割原料钢管的总根数最少的目标需要修正如下：

以上目标函数可以理解为：所采用的切割模式下所用的原料钢管总数。

约束条件为满足客户的要求，则有：

其中，dj表示用户对于第j种类型钢管的需求数，j=1,2,3,4分别表示长度为4m、5m、6m、8m的钢管。yij表示表3-7中第i种模式、第j种类型钢管的根数。

下面介绍一种可带有普遍性、可以同时确定切割模式和切割计划的方法。同问题（1）类似，一个合理的切割模式的余料不应该大于或等于客户需要的钢管的最小尺寸（本题为4m）。切割计划中只使用合理的切割模式，而由于本题中参数都是整数，所以合理的切割模式的余量不能大于3m。此外，这里仅选择总根数最少为目标进行求解。

决策变量 由于不同切割模式不能超过3种，可以用xi表示按照第i 种模式（i=1,2, 3）切割的原料钢管的根数，显然它们应当是非负整数。设所使用的第i种切割模式下每根原料钢管生产长度为4m、5 m、6 m和8m的钢管数量分别为r1i, r2i, r3i, r4i（非负整数）。

决策目标 以切割原料钢管的总根数最少为目标，即目标为：

min x1+x2+x3

约束条件 为满足客户的需求，应有

每一种切割模式必须可行、合理，所以每根原料钢管的成品量不能超过19m，也不能小于16m（余料不能大于3m），于是：

与LINDO相比，LINGO软件主要具有两大优点：除具有LINDO的全部功能外，还可用于求解非线性规划问题，包括非线性整数规划问题；内置建模语言，允许以简练、直观的方式描述较大规模的优化问题，所需的数据可以以一定格式保存在独立的文件中。

LINGO不询问是否进行敏感性分析，敏感性分析需要将来通过修改系统选项启动敏感性分析后，再调用“REPORT|RANGE”菜单命令来实现。现在同样可以把模型和结果报告保存在文件中。

一般来说，LINGO中建立的优化模型可以由五个部分组成，或称为五“段”（SECTION）：

1．集合段（SETS）：以“SETS”开始，“ENDSETS”结束，定义必要的集合变量（SET）及其元素（MEMBER，含义类似于数组的下标）和属性（ATTRIBUTE，含义类似于数组）。

2．目标与约束段：目标函数、约束条件等，没有“段”的开始和结束标记，因此实际上就是除其他四个段（都有明确的段标记）外的LINGO模型。这里一般要用到LINGO的内部函数，尤其是与集合相关的求和函数@SUM和循环函数@FOR等。

3．数据段（DATA）：以“DATA”开始，“ENDDATA”结束，对集合的属性（数组）输入必要的常数数据。格式为：“attribute（属性）=value_list（常数列表）; ”常数列表（value_list）中数据之间可以用逗号“, ”分开，也可以用空格分开（回车等价于一个空格）。

4．初始段（INIT）：以“INIT”开始，“ENDINIT”结束，对集合的属性（数组）定义初值（因为求解算法一般是迭代算法，所以用户如果能给出一个比较好的迭代初值，对提高算法的计算效果是有益的）。如果有一个接近最优解的初值，对LINGO求解模型是有帮助的。定义初值的格式为：“attribute（属性）=value_list（常数列表）; ”

5．计算段（CALC）：以“CALC”开始，“ENDCALC”结束，对一些原始数据进行计算处理。在实际问题中，输入的数据通常是原始数据，不一定能在模型中直接使用，可以在这个段对这些原始数据进行一定的“预处理”，得到模型中真正需要的数据。

普通LINGO程序设计：

    min=x1+x2+x3；
    x1*r1 1+x2*r1 2+x3*r1 3＞=5 0；
    x1*r2 1+x2*r22+x3*r23＞=1 0；
    x1*r3 1+x2*r32+x3*r33＞=2 0；
    x1*r41+x2*r42+x3*r43＞=15；
    4*r11+5*r21+6*r31+8*r41＜=19；
    4*r12+5*r22+6*r32+8*r42＜=19；
    4*r13+5*r23+6*r33+8*r43＜=19；
    4*r1 1+5*r2 1+6*r3 1+8*r4 1＞=1 6；
    4*r12+5*r22+6*r32+8*r42＞=16；
    4*r13+5*r23+6*r33+8*r43＞=16；
    x1+x2+x3＞=26；
    x1+x2+x3＜=3 1；
    @gin（x1）; @gin（x2）; @gin（x3）; @gin（r11）; @gin（r12）; @gin（r13）; @gin（r21）; @gin
（r22）; @gin（r23）; @gin（r31）; @gin（r32）; @gin（r33）; @gin（r41）; @gin（r42）; @gin
（r43）；
    End

运行结果一：

    Local optimal solution found.
    Objective value:                        28.00000
    Extended solver steps:                       29 1
    Total solver iterations:                     9681

建模化语言程序设计：

    Model：
     sets：
    number/1..3/:x；
    modes/1..4/:a, b；
    mode（number, modes）:r；

    endsets
    data：
    a=50,10,20,15；
    b=4,5,6,8；
    enddata
    min=@sum（number∶x）；
    @for（modes（i）:@sum（number（j）:x（j）*r（j, i））＞=a（i））；
    @for（number（i）:@sum（modes（j）:b（j）*r（i, j））＜=19）；
    @for（number（i）:@sum（modes（j）:b（j）*r（i, j））＞=1 6）；
    @sum（number∶x）＞=26；
    @sum（number∶x）＜=3 1；
    @for（number:@gin（x））；
    @for（mode:@gin（r））；
    End

运行结果二：

    Local optimal solution found.
    Objective value:                        28.00000
    Objective bound:                        28.00000
    Infeasibilities:                        0.000000
    Extended solver steps:                   166
    Total solver iterations:                  11142

例3.4 太阳影子定位问题

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。

1．建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00—15:00之间天安门广场（北纬39°54′26″，东经116°23′29″）3米高的直杆的太阳影子长度的变化曲线。

2．根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于表3-8的影子顶点坐标数据，给出若干个可能的地点。

问题分析

问题一：本小题要求建立影子长度变化模型，并将该模型应用于实际地理背景，进行太阳影子的仿真计算。考虑引起影子长度变化的直接因素是太阳高度角和杆子测量高度，据此分析影响太阳高度角变化的因子，最终得到影响影长变化的因子有太阳高度角、太阳赤纬角、太阳方位角和时角等参数。将模型应用到题目所给时间和位置上进行求解，绘制不同时刻下直杆太阳影子长度的变化曲线。

问题二：本小题要求根据已知某日期和时刻所对应的位置坐标，建立模型并应用到实际数据中，找出满足要求的地理位置。考虑所给数据中坐标系朝向未知，分析造成坐标变化的因素，考虑轨迹随时间的变化，基于坐标变化思想，对任意时刻的方位角求解。用最小二乘法使真实轨迹和计算轨迹偏差最小。

模型设计

问题一

考虑任意物体影子长度都是由太阳光照射引起的，且物体高度和太阳高度角是影响影长的两个直接因素。某一地点太阳高度角在不同时刻是变化的，由于地球存在公转，且地球地心一直在赤道平面上，故太阳直射点在南北回归线之间移动，进而要考虑物体所在地理纬度。考虑地球绕太阳公转时，地轴指向不变，导致太阳和地心连线与赤道平面夹角发生变化，故需考虑太阳赤纬角和太阳方位角。考虑地球自身在自转，故在探讨影子长度几何关系时还需要考虑时角等因素。

综上分析，在研究影子长度相关的几何参数时，考虑对太阳方位角、太阳赤纬角、太阳高度角和时角等参数进行分析定位，各参数的几何意义分别为：

太阳赤纬角δ：太阳中心和地球中心的连线与地球赤道平面的夹角。其中，在春秋分时刻夹角最小为0°，在夏至和冬至时角度达到最大为±23°26′36″；

太阳高度角h：地球表面上任意一点和太阳的连线与地平线的夹角；

时角t：单位时间内地球自转的角度，定义正午时角为0°，记上午时角为负值，下午时角为正值。

太阳赤纬角是地球赤道面与日地中心连线的夹角，黄赤交角为黄道面与赤道面的交角，是赤纬角的最大值。查阅相关文献，目前地球的黄赤交角约为23°26′，在日地二体系统中，选定地球为参照物，假设太阳沿赤道面绕地球做匀速圆周运动。定义积日零点为2015年1月1日，以2015年春分日为3月21日，与2015年1月1日相差79天。由于春分时δ=0°，因此以春分日为基准，则得到赤纬角的计算式为：

考虑太阳的运行轨迹，计算太阳对地球上某一物体的影子长度变化轨迹，查阅文献得到主要由当地的地理纬度和时间这两个因素决定。假设物体所在地理纬度记为φ，地理经度记为α，北京时间记为t0，北京时间所在经度记为α0，各参数的意义如下：

时角t：考虑任一经度位置的物体在计算时角时都是以北京时间为基准，且规定正午时角为0°。故得到时角计算式为：

太阳高度角h：入射至地表某点的太阳光线与该点切平面夹角，计算公式为：

h=sin-1（sinφsinδ+cosφcosδcost）, h∈[-90°,90°]

太阳赤纬角δ：太阳中心和地球中心的连线与地球赤道平面的夹角，由赤纬角的理论推导得到公式为：

综上分析得到，影子长度函数式为：

l=H×coth=H×cot（sin-1（sinφsinδ+cosφcosδcost））

考虑在影长模型中根据经纬度和时间对各个时刻的太阳高度角求解，进而求得影子长度。现将该模型应用到实际地理背景中，已知拍摄地点为天安门广场（北纬39°54′26″，东经116°23′29″）时，且拍摄时间为2015年10月22日北京时间9:00—15:00。此时得到杆长为3m的物体在各整点时的影子长度如表3-9所示。

分析表3-9中数据可得，在该时段内呈现的影子长度先减小后增大，且总大于杆子实际长度。绘制影子长度随时间连续变化曲线如图3-6所示。

分析图3-6中影子长度变化轨迹得到，影子长度在北京时间12:14时刻影子长度最短为3.7992m，即此时刻太阳高度角最大为38.296°。影子长度变化率先减小后增大，即越靠近北京时间12:14时刻，太阳高度角和影子长度单位时间的变化率越小，且曲线是连续变化的。联系实际背景可得，在北京时间10月22日，太阳直射点在南半球，且往南回归线方向移动，故相对应天安门这一地点的影子长度均大于实际长度是合理的，一定程度上也说明了模型的可靠性。

问题二

考虑对给定太阳顶点影子坐标数据，建立模型确定所观测点的地理位置。建立确定直杆所处位置与时间的数学模型，由于描述太阳状态的角度对时间十分敏感，无法直接使用北京时间，为了计算直杆所处地点，需要转换为真太阳时，即当地的地方时间，建立优化模型，运用最小二乘法使得预测影长轨迹与实际影长轨迹之间误差最小，进而得到可能的观测位置。

由问题一模型几何求解得到影长与太阳高度角之间函数关系为：

l=H×coth

其中，h为太阳高度角，h=arcsin（sinφsinδ+cosφcosδcost）。

考虑建立平面直角坐标系xOy, x轴正方向为正东，y轴正方向为正北，则影长在x轴和y 轴的分量分别为：

其中，A为太阳方位角。

考虑表3-8中所给出的x和y 坐标数据没有规定坐标轴方向，由于任意两个时刻的坐标系角度是变化的，故需要统一坐标系。再用最小二乘法计算求解，将得到的影长坐标数据与真实坐标数据进行比较。坐标变换示意图如图3-7所示。

将各时刻坐标系旋转β角后，由图3-7转换示意图分析得到新旧坐标变换公式为：

其中，β为两平面直角坐标系变化前后旋转角度；x和y 为变化前坐标；x′和y′为变化后坐标。

求解的目标是得到若干个可能的直杆所在地点，为使得阴影的长度和位置即轨迹与附件所给尽可能相近，利用最小二乘法建立如下优化模型。目标函数为：

其中，xi和yi为真实的影子横纵坐标；x′i和y′i为预测得到的影子横纵坐标。

约束条件为：

参数关系为：

对上述几何关系式中未知参数在全范围遍历使得计算得到轨迹与真实轨迹之间相对误差尽可能小，并运用最小二乘法，通过MATLAB对优化模型进行计算求解。通过对优化模型求解，得到满足影子轨迹最优和影子长度最优目标下的位置信息如表3-10所示。考虑优化模型求解后，得到通过影长优化的地理位置位于海南省。