【数学建模算法】（6）非线性规划：定义和实例

前面的几节基本都围绕着线性规划和可转化为线性规划的问题来介绍，这一节开始我们将介绍非线性规划

1.非线性规划

1.1.实例和定义

如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。

下面通过实例归纳出非线性规划数学模型的一般形式，介绍有关非线性规划的基本概念。

例1 （投资决策问题）某企业有 $n$ 个项目可供选择，并且至少要对其中一个项目投资。已知该公司有总资金 $A$ 元，投资于第 $i$ 个项目需花费资金 $a_{i}$ 元，并预计可收益 $b_{i}$ 元。试确定最佳投资方案。

解：设投资决策变量为：
$x_{i}=\left\{\begin{array}{c}{1} ，决定投资第i个项目\\ {0}，决定不投资第i个项目\end{array}\right.$
则投资总额为 $\sum_{i=1}^{n} a_{i} x_{i}$ ，投资总收益为 $\sum_{i=1}^{n} b_{k i}$ ，考虑到限制条件：至少要对一个项目投资，并
且总的投资金额不能超过总资金 $A$ ，所以有限制条件如下：
$0<\sum_{i=1}^{n} a_{i} x_{i} \leq A$
因为 $x_{i}$ 只能取0或1，所以有非线性限制条件：
$x_{i}\left(1-x_{i}\right)=0, i=1, \cdots, n$
本个问题对于目标函数的确定我们可以这样考虑：投资决策问题的目标是让成本最小的同时利润最大，所以我们极大化总收益与总投资之比。

所以目标函数为：
$\max Q=\frac{\sum_{i=1}^{n} b_{i} x_{i}}{\sum_{i=1}^{n} a_{i} x_{i}}$
限制条件为：
s. t. $\quad 0<\sum_{i=1}^{n} a_{i} x_{i} \leq A$
$x_{i}\left(1-x_{i}\right)=0, i=1, \cdots, n$

这就是上面这个例题的完整模型，上面这个例题中，目标函数和限制条件中都有非线性函数，当目标函数和限制条件中至少有一个非线性时，这类问题称为非线性规划问题，可概括为一般形式。

$\min f(x)$
$\begin{array}{ll}{\text { S.t. }} & {h_{j}(x) \leq 0, \quad j=1, \cdots, q} \\ {} & {g_{i}(x)=0, \quad i=1, \cdots, p}\end{array}$
以上模型简称为(NP)模型，其中 $x=\left[\begin{array}{lll}{x_{1}} & {\cdots} & {x_{n}}\end{array}\right]^{T}$ 称为模型的决策变量， $f$ 称为目标函数， $g_{i}(i=1, \cdots, p)$ 和 $h_{j}(j=1, \cdots, q)$ 称为约束函数。另外， $g_{i}(x)=0 \quad(i=1, \cdots, p)$ 称为等式约束， $h_{j}(x) \leq 0 \quad(j=1, \cdots, q)$ 称为不等式约束。

1.2.线性规划和非线性规划的区别

如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。

1.3实例：飞行管理问题

在约10000m高空的某边长160km的正方形区域内，常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均又计算机记录其数据，以便进行飞行管理。当一家欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机相撞。如果会相撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，为避免碰撞。先假定条件如下：

1.不碰撞的标准位任意两架飞机的距离大于8km;
2.飞机飞行方向角调整的幅度不应超过30度
3.所有飞机飞行速度均为每小时800km
4.进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60km以上
5.最多考虑六架飞机
6.不必考虑飞机离开此区域后的情况

请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过 0.01 度），要求飞机飞行方向角调整的幅度尽量小。
设该区域 4 个顶点的座标为(0,0)，(160,0)，(160,160)，(0,160)。
记录数据见下表。

飞机编号	横坐标x	纵坐标y	方向角（度）
1	150	140	243
2	85	85	236
3	150	155	220.5
4	145	50	159
5	130	150	230
新进入	0	0	52

注：方向角指飞行方向与 $x$ 轴正向的夹角。
为方便之后的讨论，我们引进如下记号。
$D$ 为飞行管理区域的边长。
$\Omega$ 为飞行管理区域，取直角坐标系使其为 $[0, D] \times[0, D]$
$a$ 为飞机飞行速度， $a=800 \mathrm{km} / \mathrm{h}$
$\left(x_{i}^{0}, y_{i}^{0}\right)$ 为第 $i$ 架飞机的初始位置。
$\left(x_{i}(t), y_{i}(t)\right)$ 为第 $i$ 架飞机在 $t$ 时刻的位置
$\theta_{i}^{0}$ 为第 $i$ 架飞机的原飞行方向角，即飞行方向与 $x$ 轴夹角： $0 \leq \theta_{i}^{0}<2 \pi$
$\Delta \theta$ 为第 $i$ 架飞机的方向角调整， $-\frac{\pi}{6} \leq \Delta \theta_{i} \leq \frac{\pi}{6}$
$\theta_{i}=\theta_{i}^{0}+\Delta \theta_{i}$ 为第 $i$ 架飞机调整后的飞行航向角

1.3.1.模型构建过程

根据相对运动的观点在考察两架飞机 $i$ 和 $j$ 的飞行时，可以将飞机 $i$ 视为不动而飞机 $j$ 以相对速度
$v_{i j}=v_{j}-v_{i}=\left(a \cos \theta_{j}-a \cos \theta_{i}, a \sin \theta_{j}-a \sin \theta_{i}\right)$
相对于飞机 $i$ 运动，对上式进行和差化积得：
$\begin{aligned} v &=2 a \sin \frac{\theta_{j}-\theta_{i}}{2}\left(-\sin \frac{\theta_{j}+\theta_{i}}{2}, \cos \frac{\theta_{j}+\theta_{i}}{2}\right) \\ &=2 a \sin \frac{\theta_{j}-\theta_{i}}{2}\left(\cos \left(\frac{\pi}{2}+\frac{\theta_{j}+\theta_{i}}{2}\right), \sin \left(\frac{\pi}{2}+\frac{\theta_{j}+\theta_{i}}{2}\right)\right) \end{aligned}$
不妨设 $\theta_{j} \geq \theta_{i}$ ，此时相对飞行方向角为 $\beta_{i j}=\frac{\pi}{2}+\frac{\theta_{i}+\theta_{j}}{2}$

由于两架飞机的初始距离是：
$r_{i j}(0)=\sqrt{\left(x_{i}^{0}-x_{j}^{0}\right)^{2}+\left(y_{i}^{0}-y_{j}^{0}\right)^{2}}$

只要相对速度方向不在扇形区域内，就不会有危险

只要相对飞行方向角满足：

时，两飞机就不可能相撞。

假设 $\beta_{i j}^{0}$ 为调整前 $j$ 架飞机相对于第 $i$ 架飞机的相对速度（矢量）与这两架飞机连线（从 $j$ 指向 $i$ 的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负）。得两架飞机不碰撞的条件为：
$\left|\beta_{i j}^{0}+\frac{1}{2}\left(\Delta \theta_{i}+\Delta \theta_{j}\right)\right|>\alpha_{i j}^{0}$
其中：
$\beta_{m n}^{0}=相对速度v_{mn}的幅角-从n指向m的连线矢量的幅角$
$=\arg \frac{e^{i \theta_{n}}-e^{i \theta_{n}}}{\left(x_{m}+i y_{m}\right)-\left(x_{n}+i y_{n}\right)}$
（注意 $\beta_{m n}^{0}$ 表达式中的 $i$ 表示虚数单位）这里我们利用复数的幅角，可以很方便地计算角度 $\beta_{m n}^{0}(m, n=1,2, \cdots, 6)$

本问题中的优化目标函数可以有不同的形式：如使所有飞机的最大调整量最小；所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和最小为目标函数，可以得到如下的数学规划模型：
$\min \sum_{i=1}^{6}\left|\Delta \theta_{i}\right|$
$\begin{array}{rcrc}{\text { s.t. }} & {\left|\beta_{i j}^{0}+\frac{1}{2}\left(\Delta \theta_{i}+\Delta \theta_{j}\right)\right|>\alpha_{i j}^{0},} & {i, j=1,2, \cdots, 6, i \neq j} \\ {} & {\left|\Delta \theta_{i}\right| \leq 30^{\circ},} & {i=1,2, \cdots, 6}\end{array}$

利用matlab程序计算 $\alpha_{i j}^{0}$ 和 $\beta_{i j}^{0}$ ：

clc,clear
x0=[150 85 150 145 130 0];
y0=[140 85 155 50 150 0];
q=[243 236 220.5 159 230 52];
xy0=[x0; y0];
d0=dist(xy0); %求矩阵各个列向量之间的距离
d0(find(d0==0))=inf;
a0=asind(8./d0) %以度为单位的反函数
xy1=x0+i*y0
xy2=exp(i*q*pi/180)
for m=1:6
    for n=1:6
        if n~=m
            b0(m,n)=angle((xy2(n)-xy2(m))/(xy1(m)-xy1(n)));
        end
    end
end
b0=b0*180/pi;
dlmwrite('txt1.txt',a0,'delimiter', '\t','newline','PC');
fid=fopen('txt1.txt','a');
fwrite(fid,'~','char'); %往纯文本文件中写 LINGO 数据的分割符
dlmwrite('txt1.txt',b0,'delimiter', '\t','newline','PC','-append','roffset', 1)

matlab会输出一个txt文件，将这个txt文件输入Lingo软件。

model:
sets:
plane/1..6/:delta;
link(plane,plane):alpha,beta;
endsets
data:
alpha=@file('txt1.txt'); !需要在alpha的数据后面加上分隔符"~";
beta=@file('txt1.txt');
enddata
min=@sum(plane:@abs(delta));
@for(plane:@bnd(-30,delta,30));
@for(link(i,j)|i#ne#j:@abs(beta(i,j)+0.5*delta(i)+0.5*delta(j))>a
lpha(i,j));
end

求解的最优解 $\Delta \theta_{3}=2.83858^{\circ}, \Delta \theta_{6}=0.7908^{\circ}$ ，其他调整角度为0.

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 203,362评论 5赞 477
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,330评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,247评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,560评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,580评论 5赞 365
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,569评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,929评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,587评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,840评论 1赞 297
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,596评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,678评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,366评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,945评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,929评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,165评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 43,271评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,403评论 2赞 342