东南大学数学建模试题(含答案)

数学建模竞赛参考答案数学建模竞赛参考答案数学建模竞赛是一项旨在培养学生综合运用数学知识和解决实际问题能力的竞赛活动。

参赛者需要通过分析问题、建立数学模型、求解问题等环节，最终给出合理的答案和解决方案。

在这篇文章中，我们将为大家提供一些数学建模竞赛的参考答案，希望能够给参赛者们提供一些启示和帮助。

第一题：某公司的销售额预测问题描述：某公司希望通过过去几年的销售数据，预测未来一年的销售额。

请根据给定的销售数据，建立合适的数学模型，并给出未来一年的销售额预测值。

解答思路：根据问题描述，我们可以将销售额看作是时间的函数，即销售额随时间变化。

可以使用回归分析的方法来建立数学模型。

首先，我们将销售额作为因变量，时间作为自变量，通过拟合曲线来预测未来一年的销售额。

我们可以选择多项式回归模型来拟合曲线。

通过将时间作为自变量，销售额作为因变量，进行多项式回归分析，可以得到一个多项式函数，该函数可以描述销售额随时间变化的趋势。

然后，我们可以使用该多项式函数来预测未来一年的销售额。

将未来一年的时间代入多项式函数中，即可得到未来一年的销售额预测值。

第二题：城市交通流量优化问题描述：某城市的交通流量问题日益突出，如何优化交通流量成为了当地政府亟待解决的难题。

请根据给定的交通数据和道路拓扑结构，建立合适的数学模型，并给出交通流量优化的方案。

解答思路：根据问题描述，我们可以将城市的交通流量看作是网络中的流量分配问题。

可以使用网络流模型来建立数学模型。

首先，我们需要将城市的道路网络抽象成一个有向图，节点表示交叉口，边表示道路，边上的权值表示道路的容量。

然后，我们可以使用最小费用最大流算法来求解交通流量优化的方案。

该算法可以通过调整道路上的流量分配，使得整个网络中的流量达到最大，同时满足道路容量的限制。

通过计算最小费用最大流，可以得到交通流量优化的方案。

最后，我们可以根据最小费用最大流算法的结果，对交通流量进行合理调控。

例如，可以调整信号灯的时长，优化交通信号控制系统，减少交通拥堵现象，提高交通效率。

09级数模试题1. 把四只脚的连线呈长方形的椅子往不平的地面上一放，通常只有三只脚着地，放不稳，然后稍微挪动几次，就可以使四只脚同时着地，放稳了。

试作合理的假设并建立数学模型说明这个现象。

（15分）解：对于此题，如果不用任何假设很难证明，结果很可能是否定的。

因此对这个问题我们假设 ：（1）地面为连续曲面（2）长方形桌的四条腿长度相同（3）相对于地面的弯曲程度而言，方桌的腿是足够长的（4）方桌的腿只要有一点接触地面就算着地。

那么，总可以让桌子的三条腿是同时接触到地面。

现在，我们来证明：如果上述假设条件成立，那么答案是肯定的。

以长方桌的中心为坐标原点作直角坐标系如图所示，方桌的四条腿分别在A 、B 、C 、D 处，A 、B,C 、D的初始位置在与x 轴平行，再假设有一条在x 轴上的线ab,则ab 也与A 、B ，C 、D 平行。

当方桌绕中心0旋转时，对角线 ab 与x 轴的夹角记为θ。

容易看出，当四条腿尚未全部着地时，腿到地面的距离是不确定的。

为消除这一不确定性，令 ()f θ为A 、B 离地距离之和，()g θ为C 、D 离地距离之和，它们的值由θ唯一确定。

由假设（1），()f θ,()g θ均为θ的连续函数。

又由假设（3），三条腿总能同时着地， 故()f θ()g θ=0必成立（∀θ）。

不妨设(0)0f =,(0)0g >g （若(0)g 也为0，则初始时刻已四条腿着地，不必再旋转），于是问题归结为：已知()f θ,()g θ均为θ的连续函数，(0)0f =,(0)0g >且对任意θ有00()()0f g θθ=，求证存在某一0θ，使00()()0f g θθ=。

证明：当θ=π时，AB 与CD 互换位置，故()0f π>，()0g π=。

作()()()h f g θθθ=-，显然，()h θ也是θ的连续函数，(0)(0)(0)0h f g =-<而()()()0h f g πππ=->，由连续函数的取零值定理，存在0θ，00θπ<<，使得0()0h θ=，即00()()f g θθ=。

东 南 大 学 考 试 卷（A 卷）课程名称 数学建模与数学实验 考试学期 2010-2011-2 得分 适用专业 各专业 考试形式 闭卷 考试时间长度 120分钟 （考试可带计算器） 所有数值结果精度要求为保留小数点后两位 一．填空题：（每题2分，共10分） 1. 用Matlab 做AHP 数学实验，常用的命令有 ， 等等。

2. 矩阵A 关于模36可逆的充要条件是： 。

3. 泛函332230()()2()3J x x t t x t t dt ⎡⎤=++⎣⎦⎰&取极值的必要条件为 。

4. 请补充一致矩阵缺失的元素136A ⎛⎫ ⎪= ⎪ ⎪⎝⎭。

5. 请列出本人提交的上机实验内容（标题即可） 。

二．选择题：（每题2分，共10分） 1. 在下列Leslie 矩阵中，能保证主特征值唯一的是 ( ) A. 0230.20000.40⎛⎫ ⎪ ⎪ ⎪⎝⎭； B. 0 1.200.10000.30⎛⎫ ⎪ ⎪ ⎪⎝⎭； C. 0070.30000.10⎛⎫ ⎪ ⎪ ⎪⎝⎭； D.以上都对 2. 下列论述正确的是 （ ） A.判断矩阵一定是一致矩阵 B.正互反矩阵一定是判断矩阵 C.能通过一致性检验的矩阵是一致矩阵 D.一致矩阵一定能通过一致性检验 3. n 阶Leslie 矩阵有 个零元素。

（ ）A.不超过2(1)n -；B.不少于2(1)n -；C.恰好2(1)n -；D.恰好21n -4. Matlab 软件内置命令不可以 （ ）A.求矩阵的主特征值B. 做曲线拟合;C. 求解整数线性规划D. 求样条插值函数5. 关于等周问题，下面的描述不正确的有 （ ）A.目标泛函可以表示为最简泛函；B.条件泛函为最简泛函；C.条件泛函取值为常数；D. 函数在区间两个端点处可以取任意值三．判断题（每题2分，共10分）1. 马氏链模型中，矩阵一定有特征值1。

（ ）2. 插值函数不要求通过样本数据点。

《专题数学建模与数学探究》试卷(答案在后面)一、单选题（本大题有8小题，每小题5分，共40分）1、答案：B解析：在数学建模中，对于给定的实际问题，通过建立数学模型来解决问题的方法属于哪种类型？A. 线性规划B. 非线性规划C. 概率统计D. 图论2、答案：C2、某公司生产某种产品，已知生产该产品的成本函数为(C(x)=50x+1000)（单位：元），其中(x)表示生产数量（单位：个）。

若该公司希望将每件产品的售价定为成本加成30%，则为了确保公司的总利润不低于1500元，至少需要生产多少件产品？A. 50B. 60C. 70D. 803、某地气温随时间变化的规律可以用函数(T(t)=−0.2t2+4t+15)来描述，其中(T)表示温度（℃），(t)表示时间（小时）。

若从上午9点开始测量，则到中午12点时，该地的平均温度是多少？A. 21.5℃B. 22℃C. 22.5℃4、某公司为了研究某种产品的销售情况，收集了过去一年的月销售额数据，并通过线性回归模型预测未来几个月的销售额。

若已知过去12个月的月销售额分别为10万元、12万元、15万元、18万元、20万元、23万元、26万元、29万元、32万元、35万元、38万元、40万元。

利用这些数据训练得到的线性回归方程为(y=1.5x+5)，其中(x)代表月份，(y)代表销售额（单位：万元）。

请问，根据该模型预测，第13个月的销售额预计是多少？A. 41万元 B) 42万元 C) 43万元 D) 44万元5、某商店为了预测某种商品的需求量，收集了过去几年该商品的销售数据，发现其销售量y与时间x（年份）之间存在线性关系，经过分析得出回归方程为(y=100+5x)。

若预测2023年的销售量，则预计销售量为：A. 600 B) 700 C) 800 D) 9006、某工厂生产一种产品，其成本函数为(C(x)=200+5x)元，其中(x)是生产的数量。

若该工厂希望以每件产品(30)元的价格销售，为了保证不亏本，至少需要生产多少件产品？A. 40 件B. 50 件C. 60 件D. 70 件7、一辆汽车在平直公路上匀速行驶，若该汽车的功率为P，行驶速度为v，则根据公式P=mvv（m为汽车质量，v为速度），可以计算出汽车所受阻力f的大小为（）。

实用标准文案精彩文档东南大学考试卷（A卷）姓名学号班级课程名称数学建模与实验考试学期 09-10-2 得分适用专业各专业考试形式闭卷考试时间长度120分钟实用标准文案精彩文档一．填空题：（每题2分，共10分）1. 阻滞增长模型0.5(10.001)(0)100dx x x dtx ⎧=-⎪⎨⎪=⎩的解为 。

2. 用Matlab 做常微分方程数学实验，常用的命令有 。

3. 整数m 关于模12可逆的充要条件是： 。

4. 根据Malthus 模型,如果自然增长率为2%，则人口数量增长为初值3倍所需时间为(假设初值为正) 。

5. 请补充判断矩阵缺失的元素13192A ⎛⎫⎪=⎪ ⎪⎝⎭。

二．选择题：（每题2分，共10分）1. 在下列Leslie 矩阵中，不能保证模最大特征值唯一的是 ( )A. 0230.20000.40⎛⎫ ⎪ ⎪ ⎪⎝⎭；B. 1.1 1.230.20000.40⎛⎫ ⎪ ⎪ ⎪⎝⎭；C. 0030.20000.40⎛⎫ ⎪⎪ ⎪⎝⎭； D.以上都不对2. 判断矩阵能通过一致性检验的标准是 （ )A. 0.1CR <B. 0.1CI <C. 0.1CR >D.0.01CR < 3. 模28倒数表中可能出现的数是 ( )A. 12B.5C.14D.74. 线性最小二乘法得到的函数不可能为 （ ）A.线性函数B. 对数函数C. 样条函数D. 指数函数5. 关于泛函极值问题，下面的描述正确的有 （ ） A.泛函()J x 在x *处取极值的充要条件是泛函变分()0J x δ*=； B. 泛函()J x 在x *处取极值的充分条件是泛函变分()0J x δ*=； C. 泛函()J x 在x *处取极值的必要条件是泛函变分()0J x δ*=；D. A,B,C 均正确三．判断题（每题2分，共10分） 1. Hill 密码体系中，任意一个可逆矩阵都可以作为加密矩阵。

（ ） 2. 拟合函数不要求通过样本数据点。

数学建模模拟试题及答案一、填空题(每题 5 分，共 20 分)1.一个连通图能够一笔画出的充分必要条件是.2. 设银行的年利率为 0.2，则五年后的一百万元相当于现在的万元.3. 在夏季博览会上，商人预测每天冰淇淋销量N 将和下列因素有关：(1) 参加展览会的人数n； (2)气温T 超过10o C；(3)冰淇淋的售价p .由此建立的冰淇淋销量的比例模型应为 .4. 如图一是一个邮路，邮递员从邮局 A 出发走遍所有 A长方形街路后再返回邮局 .若每个小长方形街路的边长横向均为 1km，纵向均为 2km，则他至少要走 km .二、分析判断题(每题 10 分，共 20 分)1. 有一大堆油腻的盘子和一盆热的洗涤剂水。

为尽量图一多洗干净盘子，有哪些因素应予以考虑？试至少列出四种。

2. 某种疾病每年新发生 1000 例，患者中有一半当年可治愈 .若 2000 年底时有1200 个病人，到 2005 年将会出现什么结果？有人说，无论多少年过去，患者人数只是趋向 2000 人，但不会达到 2000 人，试判断这个说法的正确性 .三、计算题(每题 20 分，共 40 分)1. 某工厂计划用两种原材料A, B 生产甲、乙两种产品，两种原材料的最高供应量依次为 22 和 20 个单位；每单位产品甲需用两种原材料依次为 1 、1 个单位，产值为 3 (百元)；乙的需要量依次为 3、1 个单位，产值为 9 (百元)；又根据市场预测，产品乙的市场需求量最多为 6 个单位，而甲、乙两种产品的需求比不超过 5： 2，试建立线性规划模型以求一个生产方案，使得总产值达到最大，并由此回答：(1) 最优生产方案是否具有可选择余地？若有请至少给出两个，否则说明理由 .(2) 原材料的利用情况 .2. 两个水厂A1 , A2将自来水供应三个小区B1 , B2 , B3 , 每天各水厂的供应量与各小区的需求量以及各水厂调运到各小区的供水单价见下表 .试安排供水方案，使总供水费最小？四、 综合应用题(本题 20 分)某水库建有 10 个泄洪闸，现在水库的水位已经超过安全线，上游河水还在不断地流入 水库.为了防洪，须调节泄洪速度 .经测算，若打开一个泄洪闸， 30 个小时水位降至安全线， 若打开两个泄洪闸， 10 个小时水位降落至安全线 .现在，抗洪指挥部要求在 3 个小时内将水 位降至安全线以下，问至少要同时打开几个闸门？试组建数学模型给予解决 .注：本题要求按照五步建模法给出全过程 .小区 单价/元水厂A1A供应量 / t170B34B11 07 1B26数学建模 06 春试题模拟试题参考解答一、填空题(每题 5 分，共 20 分)1. 奇数顶点个数是 0 或 2；2. 约 40.1876 ；3. N = Kn(T10) / p, (T > 10 0 C), K 是比例常数； 4. 42.二、分析判断题(每题 10 分，共 20 分)1. 解： 问题与盘子、水和温度等因素直接相关，故有相关因素：盘子的油腻程度，盘子的温度，盘子的尺寸大小；洗涤剂水的温度、浓度； 刷洗地点 的温度等.注：列出的因素不足四个，每缺一个扣 2.5 分。

共7页 第1页东 南 大 学 考 试 卷课程名称 数学建模与实验 考试学期 05-06-3得分适用专业 各专业考试形式开卷考试时间长度 120分钟. 雇佣关系模型（18`）考虑雇员一天的工作时间t 与工资w ,回答下列问题： （1）假设雇员的满意度曲线形式为()5223w t c =+，画出每天雇员的无差别曲线族示意图。

(2）如果雇主付计时工资，求出雇主与雇员的协议曲线。

（3）如果协议为雇员每天工作10小时，给出雇员的满意度曲线。

解：1）要求体现曲线族，曲线明显具有单调递增、下凸的特点2）设(),T W 为协议曲线上的任意一点，由于雇主的计时工资曲线与雇员的满意度曲线在该点()32()5()3W kT w T k w T Tc =='==+因此()()35225233T c T T c +=+，从而解得32c T =，所以，协议曲线为：()522532W T T ⎛⎫= ⎪⎝⎭3）因为3152c T ==，所以雇员的满意度曲线为()522153w t =+18`）甲方截获了一段密文: DXYPUFUPIKUPUAUP经分析这段密文是用Hill 2密码编译的,且经分析密文字母UPUA 依次代表明文字母ANJI 。

采26个字母与0～25之间建立一一对应关系，即A ～1，B ～2，…，Z ～0, （1）求出加密矩阵A,并在模26意义下，求矩阵A 的逆矩阵。

（2）利用（1）的解密矩阵将该段密文翻译成明文解：1） 密文字母UPUA 构成密文矩阵 2121161Q ⎡⎤=⎢⎥⎣⎦,明文字母ANJI 构成明文矩阵110149P ⎡⎤=⎢⎥⎣⎦。

由Hill 加密原理可知道，(mod 26)Q AP =，因此共7页 第2页11212111*********(mod 26)(mod 26)(mod 26)161149161141191410212191611025892517 25(mod 26)(mod 26)1611213900642503A QP ---⎡⎤⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥⎢⎥-⨯-⨯⎣⎦⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦11173********(mod 26)(mod 26)9(mod 26)030101093A ---⎡⎤⎡⎤⎡⎤⎡⎤====⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦密文: DXYPUFUPIKUPUAUP 对应明文ZHENGBANRUANJIAN 三、（18`）在河流中长时间游泳会损失大量热量，因此游泳者需要选择最佳路径。

综合题目参考答案1. 赛程安排(2002年全国大学生数学建模竞赛D 题) (1)用多种方法都能给出一个达到要求的赛程.(2)用多种方法可以证明n 支球队“各队每两场比赛最小相隔场次r 的上界”(如n =5时上界为1)是⎥⎦⎤⎢⎣⎡-23n ,如: 设赛程中某场比赛是i ,j 两队, i 队参加的下一场比赛是i ,k 两队(k ≠j ),要使各队每两场比赛最小相隔场次为r ,则上述两场比赛之间必须有除i ,j ,k 以外的2r 支球队参赛,于是32+≥r n ,注意到r 为整数即得⎥⎦⎤⎢⎣⎡-≤23n r . (3)用构造性的办法可以证明这个上界是可以达到的,即对任意的n 编排出达到该上界的赛程.如对于n =8, n =9可以得到:1A 2A 3A 4A 5A 6A 7A 8A每两场比赛相隔场次数相隔场次总数 1A × 1 5 9 13 17 21 25 3,3,3,3,3,3 18 2A 1 × 20 6 23 11 26 16 4,4,4,3,2,2 19 3A5 20 × 24 10 27 15 2 2,4,4,4,3,2 19 4A 96 24 × 28 24 3 19 2,2,4,4,4,3 19 5A 13 23 10 28 × 4 187 2,2,2,4,4,4 18 6A 17 11 27 14 4 × 8 22 3,2,2,2,4,4 17 7A 21 26 15 3 18 8 × 12 4,3,2,2,2,4 17 8A25 1621972212×4,4,3,2,2,2171A2A3A 4A5A 6A 7A 8A 9A每两场比赛相隔场次数 相隔场 次总数1A × 36 6 31 11 26 16 21 1 4,4,4,4,4,4,4, 28 2A 36 × 2 27 7 22 12 17 32 4,4,4,4,4,4,3 27 3A6 2 × 35 15 30 20 25 10 3,3,4,4,4,4,4 26 4A 31 27 35 × 3 18 8 13 23 4,4,4,4,3,3,3 25 5A 11 7 15 3 × 34 24 29 19 3,3,3,3,4,4,4 24 6A 26 22 30 18 34 × 4 9 14 4,4,3,3,3,3 23 7A 16 12 20 8 24 4 × 33 28 3,3,3,3,3,3,4 22 8A 21 17 25 13 29 9 33 × 5 3,3,3,3,3,3,3, 21 9A13210231914285×3,4,3,4,3,4,324可以看到, n =8时每两场比赛相隔场次数只有2,3,4, n =9时每两场比赛相隔场次数只有3,4,以上结果可以推广,即n 为偶数时每两场比赛相隔场次数只有22-n ,12-n ,2n,n 为奇数时只有23-n ,21-n . (4)衡量赛程优劣的其他指标如平均相隔场次 记第i 队第j 个间隔场次数为ij c ,2,2,1,,,2,1-==n j n i ,则平均相隔场次为∑∑=-=-=n i n j ij c n n r 121)2(1r 是赛程整体意义下的指标,它越大越好.可以计算n =8,n =9的r ,并讨论它是否达到上界.相隔场次的最大偏差 定义||,r c Max f ij j i -=∑-=--=21|)2(|n j ij r n c Max gf 为整个赛程相隔场次的最大偏差,g 为球队之间相隔场次的最大偏差,它们都是越小越好.可以计算n =8,n =9的f ,g ,并讨论它是否达到上界.参考文献工程数学学报第20卷第5期2003 2. 影院座位设计建立满意度函数),(βαf ,可以认为α和β无关, ()()βαβαh g f -=),(,g ,h 取尽量简单的形式,如αα=)(g ;0)(=βh (030≤β),0)(h h =β)30(0>β. (1)可030≤β将作为必要条件,以α最大为最佳座位的标准.在上图中以第1排座位为坐标原点建立坐标轴x ,可以得到⎪⎭⎫⎝⎛+----⎪⎭⎫ ⎝⎛+--=⎪⎭⎫⎝⎛+--=d x x h c H d x x c H d x x c H θθαθβtan arctan tan arctan ,tan arctan β是x 的减函数.可得x ≈1.7m,即第3(或4)排处030=β.又通过计算或分析可知α也是x 的减函数,所以第3(或4)排处是最佳座位.(2)设定一个座位间隔l (如0.5m), x 从0(或030≤β处)到d D -按l 离散,对于)20~0(00θ计算α的平均值,得020=θ时其值最大.(3)可设地板线是x 的二次曲线2bx ax +,寻求a ,b 使α的平均值最大. 实际上,还应考虑前排不应挡住后排的视线.3.节水洗衣机(1996年全国大学生数学建模竞赛B 题)该问题不要求对洗衣机的微观机制(物理、化学方面)深入研究,只需要从宏观层次去把握.宏观上洗衣的基本原理是用洗涤剂通过漂洗把吸附在衣物上的污物溶于水中,再脱去污水带走污物;洗衣的过程是通过“加水——漂洗——脱水”程序的反复运行,使残留在衣物的污物越来越少,直到满意的程度;洗涤剂也是不希望留在衣物上的东西,可将“污物”定义为衣物上原有污物与洗涤剂的总和.假设每轮漂洗后污物均匀地溶于水中;每轮脱水后衣物含水量为常数c .0x ~初始污水量,~k u 第k 轮加水量,k x ~第k 轮脱水量),,2,1( =k .设每轮脱水前后污物在水中的浓度不变.于是cx c u x c xc u x c x u x n n n =+=+=--11221110,,, , 得到)()(210c u c u u c x x n n n ++=. 在最终污物量与初始污物量之比0/x x n 小于给定的清洁度条件下,求各轮加水量k u ),,1(n k =,使总用水量最小,即∑=nk k u u Min k 1()ε<++)(..21c u c u u c t s n n等价于)()(21c u c u u Min n u k +++++α=++)()(..21c u c u u t s na 为常数可得c u c u u n +==+= 21,即第n ~2轮加水量u u k =(常数),第1轮加水量c u u +=1.令cx u =,问题简化为nx Min u n ,ε<⎪⎭⎫ ⎝⎛+nx t s 11.. 其解为0→x ,即0→u ,而∞→n .这与实际上是不合理的.应该加上对u 的限制:21v u v ≤≤.则得max min n n n ≤≤,其中 max min n n n ≤≤,1)/1ln(2min +⎥⎦⎤⎢⎣⎡+=c v n α这样,n为有限的几个数,可一一比较,具体数据计算从略.参考文献:《数学的实践与认识》第27卷第1期,19974.教师工资调整方案(1995年美国大学生数学建模竞赛B 题)题目对职称提升年限表述得不甚清楚(如未提及助理教授的提升),教龄也未区分是什么职称下工作的年限,所以应该作出一些相应的简化假设.按所给信息,工资仅取决于职称和教龄.建立新方案的一种办法是将职称折合成教龄,如定义x=教龄t+7×k (对于讲师、助理教授、副教授、教授,k 分别取值0,1,2,3),然后寻求工资函数I(x),使之满足题目的要求,如I(0)=27000,I(7)=32000等,以及x 较大时022<dxId .另一种办法是职称、教龄分别对待,工资函数J(k,t)从多种函数中选择,如最简单的线性函数J(k,t)=k k k k b a t b a ,,+(k=0,1,2,3)根据一定条件确定.按照第一种办法得到的新工资方案,以职称和教龄综合指标为x 的教师的工资都应为I(x),而人们的目前工资会低于或高于它.根据题目要求,高工资不应降低,低工资则应逐渐提高,尽快达到理想值I(x).需要做的只是根据每人(目前)工资与(理想值的)差额,制定学校提供的提薪资金的分配方案.它应该是简单、合理、容易被人接受的.按以上原则可以建立不同的模型,应通过检验比较其恶劣.检验可基于题目所给数据,按照提薪计划运行若干年,考察接近理想方案的情况,即用过渡时期的情况检验模型;也可进行随机模拟,按照一定规则随机产生数据(可以包括聘用、提职、解聘、退休的人数和时间等),再按照提薪计划运行,考察接近理想方案的情况.参考文献：叶其孝，《大学生数学建模竞赛辅导教材》（四），湖南教育出版社，2001 5. 一个飞行管理问题（1995年全国大学生数学建模竞赛A 题） 设ij a 为第i 架飞机与第j 架飞机的碰撞角（即)8arcsin(ijij r a =其中ij r 为这两架飞机连线的长度），ij β为第i 架飞机相对于第j 架飞机的相对速度（矢量）与这两架飞机连线（从i 指向j 的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负），i θ为第架飞机飞行方向角调整量.本问题中的优化目标函数可以有不同的形式：如使所有飞机的最大调整量最小；所有飞机的调整量绝对值之和最小等.以所有飞机的调整量绝对值之和最小，可以得到如下的数学规划模型：∑=61i i Min θs.t. ,)(21ij j i ij a >++θθβ j i j i ≠=,6,,1,30≤i θ , 6,,1 =i为了利用LINGO 求解这个数学规划模型，可以首先采用其他数学软件计算出ij α和ij β.其实，ij α和ij β也是可以直接使用LINGO 来计算的，这相当于解关于ij α和ij β的方程，只是解方程并非LINDO 软件的特长，这里我们作为一个例子，看看如何利用LINGO 计算ij α，可输入如下模型到LINGO 求解ij α：MIDEL ： 1]SETS:2] PLANE/1..6/：x0,y0; 3] link(plane,plane):alpha,sin2: 4]ENDSETS5] @FOR(LINK(I,J)|I#NE#J:6] sin2(I,J)=64/((X0(I)-X0(J))*(X0(I)-X0(J))+ 7] (Y0(I)-Y0(J))*(Y0(I)-Y0(J))); 8] ）；9] @FOR(LINK(I,J)|I#NE#J:10] (@SIN(alpha*3./180.0))^2=SIN2; 11] ); 12]DATA:13] X0=150,85,150,145,130,0; 14] Y0=140,85,155,50,150,0; 15]endata END 计算结果如下：ija j=1 2 3 4 5 6i =1 0.000 0 5.391232.2315.091820.96342.23452 5.391 2 0.0000 4.804 0 6.61355.807 9 3.81593 32.2310 4.8040.000 0 4.364722.83372.12554 5.091 8 6.6135 4.364 7 0.0004.4.537 2.98985 20.9634 5.807922.83374.53770.000 0 2.30986 2.234 5 3.8159 2.125 5 2.98982.309 8 0.000ijβ也可类似地利用LINGO求得，计算结果如下：ijβj=1 2 3 4 5 6i =1 0.000109.263 6-128.250 024.179 8173.065 114.474 92 109.263 60.000 0-88.871 1-42.243 6-92.304 89.000 03 -128.250 0-88.871 10.00012.476 3-58.786 20.310 84 24.179 8-42.243 612.476 30.000 05.969 2-3.525.65 173.065 1-92.304 8-58.786 25.969 20.000 01.914 46 14.479.000.310 -3.5 1.910.04 9 0 0 8 256 4 4 00 0于是，该飞机管理的数学规划模型可如下输入LINGO求解：MODEL:1]SETS2] plane/1..6/:cita:3] link(plane,plane):alpha,beta;4]ENDSETS5] min=@sum(plane:@abs(cita));6] @for(plane(I):7] @bnd(-30,cita(I),30);8] );9] @fpr(link(I,j)|I#NE#J:10] @ABS(beta(I,J)+0.5*cit(I)+0.5*cita(J))11] >alpha(I,J);12] ）;13]DATA:14] A;[JA=0.000 0 5.391.2…..…2.309 8 0.000 020] ;21] BETA=0.000 010 9.263 6………1.914 4 0.000 027] ;28]enddataEND[注] alpha,beta中数据略去，见上面表格.求解结果如下：OPTIMUM FOUND AT STEP 197SOLUTION OBJECTIVE VALUE= 3.630V ARIABLE V ALUE REDUCED COSTCITA(1) 0.E-06 -1.000 000 CITA(2) -0.E-05 -0.715 033 4CITA(3) 2.557 866 1.000 000 CITA(4) -0.E-04 0.E+00 CITA(5) 0.E-05 -1.000 000 CITA(6) 1.071 594 0.E+00 ………. (以下略)由此可知最优解为：︒︒≈≈07.1,56.263θθ （其它调整角度为0）.评注：如果将目标改为最大调整量最小，则可进一步化简得到线形规划模型，也可用LINDO 或LINGO 求解.参考文献：《数学的实践与认识》第26卷第1期，1996 6. 降落伞的选择这个优化问题的决策变量是降落伞数量n 和每一个伞的半径r ，可先将n 和r 看作连续变量，建立优化模型，求得最优解后，再按题目要求作适当调整.目标函数之降落伞的费用，可以根据表1数据拟合伞面费用1C 与伞的半径r 的关系。

programi :(1) function [accum, varargout] = CircularHough_Grd(img, radrange, varargin) %Detect circular shapes in a grayscale image. Resolve their center %positions and radii.%% [accum, circen, cirrad, dbg_LMmask] = CircularHough_Grd(% img, radrange, grdthres, fltr4LM_R, multirad, fltr4accum)% Circular Hough transform based on the gradient field of an image.% NOTE: Operates on grayscale images, NOT B/W bitmaps.% NO loops in the implementation of Circular Hough transform,% which means faster operation but at the same time larger% memory consumption.%%%%%%%%% INPUT: (img, radrange, grdthres, fltr4LM_R, multirad, fltr4accum) % % img: A 2-D grayscale image (NO B/W bitmap)%% radrange: The possible minimum and maximum radii of the circles% to be searched, in the format of% [minimum radius , maximum_radius] (unit: pixels)% **NOTE**: A smaller range saves computational time and% memory.%% grdthres: (Optional, default is 10, must be non-negative)% The algorithm is based on the gradient field of the% input image. A thresholding on the gradient magnitude% is performed before the voting process of the Circular% Hough transform to remove the Uniform intensity'% (sort-of) image background from the voting process.% In other words, pixels with gradient magnitudes smaller% than 'grdthres' are NOT considered in the computation.% **NOTE**: The default parameter value is chosen for% images with a maximum intensity close to 255. For cases% with dramatically different maximum intensities, e.g.% 10-bit bitmaps in stead of the assumed 8-bit, the default% value can NOT be used. A value of 4% to 10% of the maximum% intensity may work for general cases.%% fltr4LM_R: (Optional, default is 8, minimum is 3)% The radius of the filter used in the search of local% maxima in the accumulation array. To detect circles whose% shapes are less perfect, the radius of the filter needs% to be set larger.%% multirad: (Optional, default is 0.5)% In case of concentric circles, multiple radii may be% detected corresponding to a single center position. This% argument sets the tolerance of picking up the likely% radii values. It ranges from 0.1 to 1, where 0.1% corresponds to the largest tolerance, meaning more radii % values will be detected, and 1 corresponds to the smallest % tolerance, in which case only the "principal" radius will% be picked up.%% fltr4accum: (Optional. A default filter will be used if not given)% Filter used to smooth the accumulation array. Depending % on the image and the parameter settings, the accumulation % array built has different noise level and noise pattern% (e.g. noise frequencies). The filter should be set to an% appropriately size such that ifs able to suppress the% dominant noise frequency.%%%%%%%%% OUTPUT: [accum, circen, cirrad, dbg_LMmask]%% accum: The result accumulation array from the Circular Hough% transform. The accumulation array has the same dimension % as the input image.%% circen: (Optional)% Center positions of the circles detected. Is a N-by-2% matrix with each row contains the (x, y) positions% of a circle. For concentric circles (with the same center% position), say k of them, the same center position will% appear k times in the matrix.%% cirrad: (Optional)% Estimated radii of the circles detected. Is a N-by-1% column vector with a one-to-one correspondance to the% output tircen*. A value 0 for the radius indicates a% failed detection of the circle's radius.%% dbg_LMmask: (Optional, for debugging purpose)% Mask from the search of local maxima in the accumulation % array.%%%%%%%%%% EXAMPLE #0:% rawimg = imread('TestImg_CHT_a2.bmp');% tic;% [accum, circen, cirrad] = CircularHough_Grd(rawimg, [15 60]);% toe;% figure(l); imagesc(accum); axis image;% title(,Accumulation Array from Circular Hough Transfbrm,);% figure(2); imagesc(rawimg); colormap(,gray,); axis image;% hold on;% plot(circen(:,l), circen(:,2), *r+');% for k = 1 : size(circen, 1),% DrawCircle(circen(k, 1), circen(k,2), cirrad(k), 32,，b」)；% end% hold off;% title([*Raw Image with Circles Detected% '(center positions and radii marked)*]);% figure(3); surf(accum, 'EdgeColoF, hone'); axis ij;% title('3-D View of the Accumulation Array*);%% COMMENTS ON EXAMPLE #0:% Kind of an easy case to handle. To detect circles in the image whose% radii range from 15 to 60. Default values for arguments 'grdthres',% 'fltr4LM_R', 'multirad* and ,fltr4accum, are used.%%%%%%%%%% EXAMPLE #1:% rawimg = imread('TestImg_CHT_a3.bmp');% tic;% [accum, circen, cirrad] = CircularHough_Grd(rawimg, [15 60], 10, 20);% toe;% figure(l); imagesc(accum); axis image;% title(,Accumulation Array from Circular Hough Transfbrm,);% figure(2); imagesc(rawimg); colormap('gray'); axis image;% hold on;% plot(circen(:,l), circen(:,2), T+');% for k = 1 : size(circen, 1),% DrawCircle(circen(k, 1), circen(k,2), cirrad(k), 32, 'b-');% end% hold off;% title([*Raw Image with Circles Detected% '(center positions and radii marked)*]);% figure(3); surf(accum, 'EdgeColoF, hone'); axis ij;% title(*3-D View of the Accumulation Array*);%% COMMENTS ON EXAMPLE #1:% The shapes in the raw image are not very good circles. As a result,% the profile of the peaks in the accumulation array are kind of% 'stumpy', which can be seen clearly from the 3-D view of the% accumulation array, (As a comparison, please see the sharp peaks in % the accumulation array in example #0) To extract the peak positions % nicely, a value of 20 (default is 8) is used for argument 'fltr4LM_R', % which is the radius of the filter used in the search of peaks.%%%%%%%%%% EXAMPLE #2:% rawimg = imread(,TestImg_CHT_b3 .bmp1);% fltr4img = [1 1 1 1 1; 1 2 2 2 1; 1 2 4 2 1; 1 2 2 2 1; 1 1 1 1 1];% fltr4img = fltr4img / sum(fltr4img(:));% imgfltrd = filter2( fltr4img , rawimg );% tic;% [accum, circen, cirrad] = CircularHough_Grd(imgfltrd, [15 80], 8, 10); % toe;% figure(l); imagesc(accum); axis image;% title(,Accumulation Array from Circular Hough Transfbrm,);% figure(2); imagesc(rawimg); colormap('gray'); axis image;% hold on;% plot(circen(:,l), circen(:,2), T+');% for k = 1 : size(circen, 1),% DrawCircle(circen(k, 1), circen(k,2), cirrad(k), 32, 'b-');% end% hold off;% title([*Raw Image with Circles Detected% '(center positions and radii marked)*]);%% COMMENTS ON EXAMPLE #2:% The circles in the raw image have small scale irregularities along % the edges, which could lead to an accumulation array that is bad for % local maxima detection. A 5-by-5 filter is used to smooth out the % small scale irregularities. A blurred image is actually good for the % algorithm implemented here which is based on the image's gradient % field.%%%%%%%%%% EXAMPLE #3:% rawimg = imread('TestImg_CHT_c3.bmp');% fltr4img = [1 1 1 1 1; 1 2 2 2 1; 1 2 4 2 1; 1 2 2 2 1; 1 1 1 1 1];% fltr4img = fltr4img / sum(fltr4img(:));% imgfltrd = filter2( fltr4img , rawimg );% tic;% [accum, circen, cirrad]=...% CircularHough_Grd(imgfltrd, [15 105], 8, 10, 0.7);% toe;% figure(l); imagesc(accum); axis image;% figure(2); imagesc(rawimg); colormap(,gray,); axis image;% hold on;% plot(circen(:,l), circen(:,2), *r+');% for k = 1 : size(circen, 1),% DrawCircle(circen(k, 1), circen(k,2), cirrad(k), 32,，b」)；% end% hold off;% title([*Raw Image with Circles Detected% '(center positions and radii marked)*]);%% COMMENTS ON EXAMPLE #3:% Similar to example #2, a filtering before circle detection works for% noisy image too. 'multirad* is set to 0.7 to eliminate the false% detections of the circles* radii.%%%%%%%%%% BUG REPORT:% This is a beta version. Please send your bug reports, comments and% suggestions to pengtao@ . Thanks.%%%%%%%%%%% INTERNAL PARAMETERS:% The INPUT arguments are just part of the parameters that are used by% the circle detection algorithm implemented here. Variables in the code% with a prefix ,prm_, in the name are the parameters that control the% judging criteria and the behavior of the algorithm. Default values for% these parameters can hardly work for all circumstances. Therefore, at% occasions, the values of these INTERNAL PARAMETERS (parameters that% are NOT exposed as input arguments) need to be fine-tuned to make% the circle detection work as expected.% The following example shows how changing an internal parameter could% influence the detection result.% 1. Change the value of the internal parameter 'prm LM LoBndRa* to 0.4% (default is 0.2)% 2. Run the following matlab code:% fltr4accum = [1 2 1; 2 6 2; 1 2 1];% fltr4accum = fltr4accum / sum(fltr4accum(:));% rawimg = imread(,Frame_0_0022jportion.jpg,);% tic;% [accum, circen] = CircularHough_Grd(rawimg,...% [4 14], 10, 4, 0.5, fltr4accum);% toe;% figure(l); imagesc(accum); axis image;% title(*Accumulation Array from Circular Hough Transform*);% figure(2); imagesc(rawimg); colormap(,gray,); axis image;% hold on; plot(circen(:,l), circen(:,2), "); hold off;% title('Raw Image with Circles Detected (center positions marked)*);% 3. See how different values of the parameter 'prm LM LoBndRa* could % influence the result.% Author: Tao Peng% Department of Mechanical Engineering% University of Maryland, College Park, Maryland 20742, USA% pengtao@% Version: Beta Revision: Mar. 07, 2007%%%%%%%% Arguments and parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Validation of argumentsif ndims(img)〜=2 || 〜isnumeric(img),error(*CircularHough_Grd: "img" has to be 2 dimensionaf);endif 〜all(size(img) >= 32),erro^'CircularHough Grd: "img" has to be larger than 32-by-32');endif numel(radrange)〜=2 || -isnumeric(radrange),error([*CircularHough_Grd: "radrange" has to be \ ...'a two-element vector1]);endprm_r_range = sort(max( [0,0;radrange( 1 ),radrange(2)]));% Parameters (default values)prmgrdthres = 10;prmfltrLMR = 8;prmmultirad = 0.5;funccompucen = true;funccompuradii = true;% Validation of argumentsvapgrdthres = 1;if nargin > (1 + vap_grdthres),if isnumeric(varargin{vap grdthres}) && ...varargin(vap grdthres} (1) >= 0,prm_grdthres = varargin {vapgrdthres} (1);elseerror(['CircularHough_Grd: "grdthres" has to be'a non-negative number1]);endendvap_fltr4LM = 2; % filter for the search of local maximaif nargin > (1 + vap_fltr4LM),if isnumeric(varargin{vap_fltr4LM}) && varargin{vap_fltr4LM}(1) >= 3, prmfltrLMR = varargin{vap_fltr4LM} (1);elseerror([,CircularHough_Grd: n fltr4LM_R n has to belarger than or equal to 3']);endendvap_multirad = 3;if nargin > (1 + vap multirad),if isnumeric(varargin{vap_multirad}) && ...varargin{vap multirad}(1) >= 0.1 && ...varargin {vap multirad} (1) <= 1,prmmultirad = varargin {vap_mul tirad} (1);elseerror(['CircularHough_Grd: "multirad" has to be'within the range [0.1, 1]*]);endendvap_fltr4accum = 4; % filter for smoothing the accumulation arrayif nargin > (1 + vap_fltr4accum),if isnumeric(varargin{vap_fltr4accum}) && ...ndims(varargin{vap_fltr4accum}) == 2 && ...all(size(varargin {vap_fltr4accum}) >= 3),fltr4accum = varargin {vap_fltr4accum};elseerror(['CircularHough_Grd: n fltr4accum n has to be \ ...*a 2-D matrix with a minimum size of 3-by-3']);endelse% Default filter (5-by-5)fltr4accum = ones(5,5);fltr4accum(2:4,2:4) = 2;fltr4accum(3,3) = 6;end func_compu_cen = (nargout > 1 );func_compu_radii = (nargout > 2 );% Reserved parametersdbg on = false; % debug information dbgbfigno = 4;if nargout > 3, dbg on = true; end%%%%%%%% Buildingaccumulation array %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Convert the image to single if it is not of% class float (single or double) img_is_double = isa(img, double');if ~(img_is_double || isa(img, 'single')),imgf = single(img);end% Compute the gradient and the magnitude of gradientif img_is_double,[grdx, grdy] = gradient(img);else[grdx, grdy] = gradient(imgf);endgrdmag = sqrt(grdx.A2 + grdy.A2);% Get the linear indices, as well as the subscripts, of the pixels% whose gradient magnitudes are larger than the given threshold grdmasklin = find(grdmag > prm_grdthres);[grdmask_ldxl, grdmask_IdxJ] = ind2sub(size(grdmag), grdmasklin);% Compute the linear indices (as well as the subscripts) of% all the votings to the accumulation array.% The Matlab function 'accumarray* accepts only double variable, % so all indices are forced into double at this point.% A row in matrix ,lin2accum_aJ, contains the J indices (into the % accumulation array) of all the votings that are introduced by a % same pixel in the image. Similarly with matrix linZaccum aP. rr_41inaccum = double( prm_r_range );linaccum_dr = [ (-rr_41inaccum(2) + 0.5) : -rr_41inaccum(l),... (rr_41inaccum(l) + 0.5) : rr_41inaccum(2)];lin2accum_aJ = floor(...double(grdx(grdmasklin)./grdmag(grdmasklin)) * linaccum_dr + ...repmat( double(grdmask_IdxJ)+0.5 , [ 1 ,length(linaccum_dr)])...);lin2accum_al = floor(...double(grdy(grdmasklin)./grdmag(grdmasklin)) * linaccum dr + ...repmat( double(grdmask_IdxI)+0.5 , [1 ,length(linaccum_dr)])...)；% Clip the votings that are out of the accumulation arraymask_valid_a J al =...lin2accum_aJ > 0 & lin2accum_aJ < (size(grdmag,2) + 1) & ...Iin2accum_al > 0 & lin2accum_al < (size(grdmag,l) + 1);mask_valid_aJaI_reverse =〜mask_valid_aJaI;lin2accum_aJ = lin2accum_aJ .* maskvalida J al + maskvalidaJ alreverse;lin2accum_al = lin2accum_al .* mask_valid_aJaI + mask_valid_aJaI_reverse;clear mask_valid_aJ alre verse;% Linear indices (of the votings) into the accumulation arraylin2accum = sub2ind( size(grdmag), lin2accum_al, lin2accum_aJ );lin2accum_size = size( lin2accum );lin2accum = reshape( lin2accum, [numel(lin2accum),l]);clear lin2accum_al lin2accum_aJ;% Weights of the votings, currently using the gradient maginitudes% but in fact any scheme can be used (application dependent)weight4accum =...repmat( double(grdmag(grdmasklin)) , [lin2accum_size(2), 1 ]) .* ...mask_valid_aJ al(:);clear mask_valid_aJaI;% Build the accumulation array using Matlab function 'accumarray'accum = accumarray( lin2accum , weight4accum );accum = [ accum ; zeros( numel(grdmag) - numel(accum), 1 )];accum = reshape( accum, size(grdmag));%%%%%%%% Locating local maxima in the accumulation array %%%%%%%%%%%%% Stop if no need to locate the center positions of circlesif ~func_compu_cen,return;endclear lin2accum weight4accum;% Parameters to locate the local maxima in the accumulation array% — Segmentation of 'accum' before locating LM prmuseaoi = true;prm_aoithres_s = 2;prm aoiminsize = floor(min([ min(size(accum)) * 0.25,... prm_r_range(2) * 1.5 ]));% — Filter for searching for local maxima prmfltrLMs = 1.35;prm fltrLM r = ceil( prm fltrLM R * 0.6 );prm fltrLM npix = max([ 6, ceil((prm_fltrLM_R/2)A 1.8)]);% — Lower bound of the intensity of local maximaprm LM LoBndRa = 0.2; % minimum ratio of LM to the max of'accum'% Smooth the accumulation arrayfltr4accum = fltr4accum / sum(fltr4accum(:));accum = filter2( fltr4accum, accum );% Select a number of Areas-Of^Interest from the accumulation array if prmuseaoi, % Threshold value for 'accum1prm_llm_thresl = prm_grdthres * prm_aoithres_s;% Thresholding over the accumulation array accummask = ( accum > prm llm thres 1 );% Segmentation over the mask[accumlabel, accum nRgn] = bwlabel( accummask, 8 );% Select AOIs from segmented regionsaccumAOI = ones(0,4);for k = 1 : accum nRgn,accumrgn lin = find( accumlabel = k);[accumrgn_ldxl, accumrgn_IdxJ]=...ind2sub( size(accumlabel), accumrgn lin);rgn top = min( accumrgn ldxl);rgn bottom = max( accumrgn_ldxl);rgn left = min( accumrgn ldxJ );rgn_right = max( accumrgn ldxJ );% The AOIs selected must satisfy a minimum sizeif ((rgn_right - rgn_left + 1) >= prm_aoiminsize && ...(rgn_bottom - rgn top + 1) >= prm aoiminsize ),accumAOI = [ accumAOI;...rgn top, rgn bottom, rgn left, rgn right ];endendelse% Whole accumulation array as the one AOIaccumAOI = [1, size(accum,l), 1, size(accum,2)];end% Thresholding of 'accum' by a lower boundprm LM LoBnd = max(accum(:)) * prm LM LoBndRa;% Build the filter for searching for local maxima fltr4LM = zeros(2 * prm_fltrLM_R + 1);[mesh4fLM_x, mesh4fLM_y] = meshgrid(-prm_fltrLM_R : prm fltrLM R);mesh4fLM_r = sqrt( mesh4fLM_x.A2 + mesh4fLM_y.A2 );fltr4LM_mask =...(mesh4fLM_r > prm_fltrLM_r & mesh4fLM_r <= prm fltrLM R );fltr4LM = fltr4LMfltr4LM_mask * (prm fltrLM s / sum(fltr4LM_mask(:)));if prm_fltrLM_R >= 4,fltr4LM_mask = ( mesh4fLM_r < (prm_fltrLM_r - 1))；elsefltr4LM_mask = ( mesh4fLM_r < prm fltrLM r );endfltr4LM = fltr4LM + fltr4LM mask / sum(fltr4LM_mask(:));% **** Debug code (begin)if dbg_on,dbg_LMmask = zeros(size(accum));end% **** Debug code (end)% For each of the AOIs selected, locate the local maximacircen = zeros(0,2);fbrk = 1 : size(accumAOI, 1),aoi = accumAOI(k,:); % just for referencing convenience% Thresholding of 'accum* by a lower boundaccumaoi_LBMask =...(accum(aoi(l):aoi(2), aoi(3):aoi(4)) > prm LM LoBnd );% Apply the local maxima filtercandLM = conv2( accum(aoi( 1):aoi(2), aoi(3):aoi(4)),...fltr4LM, 'same*);candLM mask = ( candLM > 0 );% Clear the margins of 'candLM mask*candLM_mask([l :prm_fltrLM_R, (end-prm_fltrLM_R+l):end], :) = 0;candLM mask(:, [l:prm_fltrLM_R, (end-prm_fltrLM_R+l):end]) = 0;% **** Debug code (begin)if dbg_on,dbg_LMmask(aoi( 1 ):aoi(2), aoi(3):aoi(4))=...dbg_LMmask(aoi( 1 ):aoi(2), aoi(3):aoi(4)) + ...accumaoi LBMask + 2 * candLM mask;end% **** Debug code (end)% Group the local maxima candidates by adjacency, compute the% centroid position for each group and take that as the center% of one circle detected[candLM label, candLM nRgn] = bwlabel( candLM_mask, 8 );fbr ilabel = 1 : candLM nRgn,% Indices (to current AOI) of the pixels in the groupcandgrp masklin = find( candLM label == ilabel);[candgrp_ldxl, candgrp_IdxJ]=...ind2sub( size(candLM label), candgrp masklin );% Indices (to 'accum') of the pixels in the groupcandgrp_ldxl = candgrp_ldxl + ( aoi(l) - 1 );candgrp IdxJ = candgrp IdxJ + ( aoi(3) - 1 );candgrp_idx2acm =...sub2ind( size(accum) , candgrp ldxl, candgrp IdxJ );% Minimum number of qulified pixels in the groupif sum(accumaoi_LBMask(candgrp_masklin)) < prm_fltrLM_npix, continue;end% Compute the centroid positioncandgrp_acmsum = sum( accum(candgrp_idx2acm));cc_x = sum( candgrp IdxJ .* accum(candgrp_idx2acm) ) / ...candgrpacmsum;cc_y = sum( candgrp_ldxl .* accum(candgrp_idx2acm) ) / ...candgrpacmsum;circen = [circen; cc_x, cc_y];endend% **** Debug code (begin)if dbg_on,figure(dbg bfigno); imagesc(dbg LMmask); axis image;title(*Generated map of local maxima1);if size(accumAOI, 1) == 1,figure(dbg_bfigno+1);surf(candLM, 'EdgeColor1, hone'); axis ij;title(,Accumulation array after local maximum filtering*);endend% **** Debug code (end)%%%%%%%% Estimation of the Radii of Circles %%%%%%%%%%%%% Stop if no need to estimate the radii of circlesif ~func_compu_radii,varargout{l} = circen;return;end% Parameters for the estimation of the radii of circlesfltr4SgnCv=[2 1 1];fltr4SgnCv = fltr4SgnCv / sum(fltr4SgnCv);% Find circle's radius using its signature curve cirrad = zeros( size(circen,l), 1 );for k = 1 : size(circen,l),% Neighborhood region of the circle for building the sgn. curve circen_round = round( circen(k,:));SCvR IO = circen_round(2) - prm_r_range(2) - 1;ifSCvR_IO<l,SCvR_I0= 1;endSCvRIl = circen_round(2) + prm_r_range(2) + 1;if SCvR Il > size(grdx,l),SCvRIl = size(grdx,l);endSCvR JO = circen round(l) - prm_r_range(2) - 1;ifSCvR_JO<l,SCvRJO = 1;endSCvRJ 1 = circenround(l) + prm_r_range(2) + 1;if SCvR Jl > size(grdx,2),SCvRJl = size(grdx,2);end% Build the sgn. curveSgnCvMat_dx = repmat( (SCvR J0:SCvR J 1) - circen(k,l),...[SCvRJl - SCvRJO +1,1]);SgnCvMat_dy = repmat( (SCvR_IO:SCvR_Il)' - circen(k,2),...[1 , SCvRJl - SCvRJO + 1]);SgnCvMat_r = sqrt( SgnCvMat dx .A2 + SgnCvMat_dy .A2 );SgnCvMatrpl = round(SgnCvMatr) + 1;f4SgnCv = abs(...double(grdx(SCvR_IO:SCvR_Il, SCvRJO:SCvRJ 1)) .* SgnCvMat_dx + ...double(grdy(SCvR_IO:SCvR Il, SCvR JO:SCvR J 1)) .* SgnCvMat dy...)./ SgnCvMat r;SgnCv = accumarray( SgnCvMat rp 1(:) , f4SgnCv(:));SgnCv_Cnt = accumarray( SgnCvMat rp 1 (:) , ones(numel(f4SgnCv), 1));SgnCv_Cnt = SgnCv_Cnt + (SgnCv_Cnt == 0);SgnCv = SgnCv ./ SgnCv_Cnt;% Suppress the undesired entries in the sgn. curve% ― Radii that correspond to short arcsSgnCv = SgnCv .* ( SgnCv_Cnt >= (pi/4 * [O:(numel(SgnCv_Cnt)-1 )]*));% ― Radii that are out of the given rangeSgnCv( 1 : (round(prm_r_range( 1))+1) ) = 0;SgnCv( (round(prm_r_range(2))+1) : end ) = 0;% Get rid of the zero radius entry in the arraySgnCv = SgnCv(2:end);% Smooth the sgn. curveSgnCv = filtfilt( fltr4SgnCv , [1] , SgnCv );% Get the maximum value in the sgn. curveSgnCv_max = max(SgnCv);if SgnCv_max <= 0,cirrad(k) = 0;continue;end% Find the local maxima in sgn. curve by 1st order derivatives% ― Mark the ascending edges in the sgn. curve as Is and% ― descending edges as OsSgnCv AscEdg = ( SgnCv(2:end) - SgnCv(l:(end-l)) ) > 0;% ― Mark the transition (ascending to descending) regionsSgnCv LMmask = [ 0; 0; SgnCv_AscEdg(l:(end-2)) ] & (〜SgnCv_AscEdg);SgnCv LMmask = SgnCvLMmask & [ SgnCv_LMmask(2:end); 0 ];% Incorporate the minimum value requirementSgnCvLMmask = SgnCvLMmask & ...(SgnCv(l:(end-l)) >= (prm_multirad * SgnCv_max));% Get the positions of the peaksSgnCv LMPos = sort( find(SgnCv_LMmask));% Save the detected radiiif isempty(SgnCvLMPos),cirrad(k) = 0;elsecirrad(k) = SgnCvLMPos(end);for i radii = (length(SgnCv LMPos) - 1) : -1 : 1,circen = [ circen; circen(k,:)];cirrad = [ cirrad; SgnCv_LMPos(i_radii)];endendend% Outputvarargout{l} = circen;varargout{2} = cirrad;if nargout > 3,varargout{3} = dbg_LMmask;endprograms：programs:2 function DrawCircle (x, y, r, nseg, S)% Draw a circle on the current figure using ploylines%% DrawCircle (x, y, r, nseg, S)% A simple function for drawing a circle on graph.%% INPUT: (x, y, r, nseg, S)% x, y: Center of the circle% r: Radius of the circle% nseg: Number of segments for the circle% S: Colors, plot symbols and line types%% OUTPUT: None%% BUG REPORT:% Please send your bug reports, comments and suggestions to% pengtao@glue. umd. edu . Thanks.% Author: Tao Peng% Department of Mechanical Engineering% University of Maryland, College Park, Maryland 20742, USA % pengtao@glue. umd. edu% Version: alpha Revision: Jan. 10, 2006theta = 0 : (2 * pi / nseg) : (2 * pi);pline_x 二r * cos(theta) + x;pline_y 二r * sin(theta) + y;plot (pline_x, pline_y, S);3function testiml二imread (' image 1. jpg');% rawimg = imread(，TestImg_CHT_c3. bmp J);rawimg=rgb2gray(iml);tic;[accum, circen, cirrad] = CircularHough_Grd(rawimg, [20 30], 5,50);circentoe;figure(1) ; imagesc(accum); axis image;title (J Accumulation Array from Circular Hough Transform，); figure (2) ; imagesc (rawimg) ; colormap (J gray，) ; axis image; hold on;plot (circen(:, 1), circen(:, 2), ' r+');for k = 1 : size (circen, 1),DrawCircle (circen (k, 1), circen (k, 2), cirrad (k), 32, ' b-'); end hold off; title(f Raw Image with Circles Detected ...'(center positions and radii marked)']);figure (3); surf(accum, ' EdgeColor，, ' none5); axis ij; title (J 3-D View of the Accumulation Array');附带图像image 1. jpg直接运行test.m即可得到上方的结果！当然方法是活的，只要合理即可行。

考试内容分布：1、线性规划2题，有1题需编程；2、非线性规划2题，有1题需编程；3、微分方程1题，需编程；4、差分方程2题，纯计算，不需编程；5、插值2题，拟合1题，纯计算，不需编程；；6、综合1题(4分)，纯计算，不需编程。

一、列出下面线性规划问题的求解模型，并给出matlab计算环境下的程序1.某车间有甲、已两台机床，可用于加工三种工件，假定这两台车床的可用台时数分别为800和900，三种工件的数量分别为400,600和500，且已知用两种不同车床加工单位数量不同工件所需的台时数和加工费用如下表。

问怎样分配车床的加工任务，才能即满足加工工件的要求，又使加工费用最低。

（答案见课本P35, 例1）2.有两个煤厂A,B，每月进煤分别不少于60t、100t，它们负责供应三个居民区的用煤任务，这三个居民区每月需用煤分别为45t, 75t, 40t。

A厂离这三个居民区分别为10km, 5km, 6km，B厂离这三个居民区分别为4km, 8km, 15km，问这两煤厂如何分配供煤，才能使总运输量最小？(1)问题分析设A煤场向这三个居民区供煤分别为x1,x2,x3;B煤场向这三个居民区供煤分别为x4,x5,x6，则min f=10*x1+5*x2+6*x3+4*x4+8*x5+15*x6，再根据题目约束条件来进行解题。

(2) 模型的求解>> f=[10 5 6 4 8 15];>> A=[-1 -1 -1 0 0 00 0 0 -1 -1 -1-1 0 0 -1 0 00 -1 0 0 -1 00 0 -1 0 0 -1];>> b=[-60;-100;-45;-75;-40];>> Aeq=[];>> beq=[];>> vlb=zeros(6,1);>> vub=[];>> [x,fval]=linprog(f,A,b,Aeq,beq,vlb,vub)Optimization terminated.(3) 结果分析x =0.0000 20.0000 40.0000 45.0000 55.0000 0.0000 fval = 960.0000即A 煤场分别向三个居民区供煤0t,20t,40t ；B 煤场分别向三个居民区供煤45t,55t,0t 可在满足条件下使得总运输量最小。