数学建模.docx - 冰豆网

资源描述

数学建模.docx

《数学建模.docx》由会员分享，可在线阅读，更多相关《数学建模.docx（33页珍藏版）》请在冰豆网上搜索。

数学建模.docx

数学建模

2002年全国大学生数学建模竞赛

（封面）

选择题号:

普通组

大专组

A√

（根据所选竞赛题目在方框内打√）

学校名称_____××××××××_____________

学生姓名×××、×××、×××

指导老师×××

全国大学生数学建模竞赛浙江赛区组委会

二〇〇二年九月

摘要：

本文主要研究了车灯线光源长度在满足光照强度的设计要求下和功率节能的最优解策略。

分别用正向光线追迹、逆向光线追迹、方程组模型求解。

得到的结果基本一致，但计算复杂度逐级下降、求解精度逐级上升，最后得出线性光源长度为4.060（mm）。

由得出的解绘出测试屏上光强分布图十分附和车灯实际照射情况。

最后，对设计规范从照射车距（安全性）、视认度（驾驶员）、车灯功率（节能原则）来评价其合理性。

关键字：

光线追迹离散化等价光源

车灯线光源的优化设计模型

×××（计算机系2000）×××（计算机系2000）×××（计算机系2000）

指导老师×××

摘要本文主要研究了车灯线光源长度在满足光照强度的设计要求和功率节能的最优解策略。

分别用正向光线追迹、逆向光线追迹、方程组模型求解。

得到的结果基本一致，但计算复杂度逐级下降、求解精度逐级上升，最后得出线性光源长度为4.060（mm）。

由得出的解绘出测试屏上光强分布图十分附和车灯实际照射情况。

最后，对设计规范从照射车距（安全性）、视认度（驾驶员）、车灯功率（节能原则）来评价其合理性。

关键字光线追迹离散化等价光源

一、问题重述

安装在汽车头部的车灯的形状为一旋转抛物面，车灯的对称轴水平地指向正前方,其开口半径36毫米，深度21.6毫米。

经过车灯的焦点，在与对称轴相垂直的水平方向，对称地放置一定长度的均匀分布的线光源。

要求在某一设计规范标准下确定线光源的长度。

该设计规范在简化后可描述如下。

在焦点F正前方25米处的A点放置一测试屏，屏与FA垂直，用以测试车灯的反射光。

在屏上过A点引出一条与地面相平行的直线，在该直线A点的同侧取B点和C点，使AC=2AB=2.6米。

要求C点的光强度不小于某一额定值（可取为1个单位），B点的光强度不小于该额定值的两倍（只须考虑一次反射）。

二、模型假设

假设1：

线光源是理想的线光源，分布均匀且无粗细。

假设2：

只考虑经过一次反射。

假设3：

所有的反射都符合反射定律。

即，入射光线与反射光线的方向与反射面的法线矢量的夹角相等，并且所有的光线都在同一个平面内。

并且属无损反射。

假设4：

电能全部转化成灯丝的光能。

假设５：

不考虑配光镜对光线传输方向上的影响

三、问题分析

由抛物面的已知数据：

开口半径36毫米，深度21.6毫米。

可以马上得出抛物面方程：

z2＋y2＝60x

（由抛物线：

y2＝60x，z＝0，和焦点坐标是（15,0），沿中轴线z轴旋转而成）用METLAB软件可以绘出抛物面的大致立体图形。

如图：

模型要求分析

根据题意要求，在测试屏上要求C点的光强度不小于某一额定值（可取为1个单位），B点的光强度不小于该额定值的两倍。

根据上图，分析抛物面对光反射的特性和由假设2可知，B点和C点的光强度取决于两种光的叠加——经过反射镜反射后射入B（或C）的光线和线光源直接射入B（或C）的光线。

问题1，从线光源发出，经过反射镜后，照射在B,C上光强度。

问题2，由线光源直接照射到测试屏上B,C点的光强度。

汽车灯内部分析

对车灯内部结构分析之后我们可以得到如下特性：

1．电压U是恒定的。

2．灯丝横截面面积为恒定极小S。

（由假设1可得）

3．在稳定发光时由于温度不变，则电阻率ρ恒定。

电阻R=ρl/S也恒定。

因为车灯发光电路是一个非纯电阻电路。

不考虑热效应，所以电源供给的能量完全转化为光幅射能。

（由假设4可得）。

由P=U2/R，我们可以推出P=U2/（ρl/S）

转换后：

Pl=U2S/ρ

∴由上1，2，3，可得U2S/ρ为常量C

由此，我们可以得出线光源与功率成反比的关系

Pl=C

根据式子要想功率尽可能的小，线光源要尽可能的大，且要满足在测试屏上B和C点的光强比要求。

四、模型的建立和求解

1，建立空间坐标系

由题目中可知，因为这个状如抛物面的车灯的开口半径为36毫米,深度为21.6。

则设车灯抛物面的对称轴为X轴，车灯的抛物面的顶点作为我们的坐标原点。

建立三维空间坐标系。

根据已知条件，设在X-Y平面里的抛物线方程为y2＝2px,且过点（21.6，36）（单位：

毫米）则，可以求得p=30，那么，抛物线为y2＝60x，z＝0。

该抛物线沿中轴线x轴旋转，可以得出车灯的抛物面：

z2＋y2＝60x。

空间图形如图：

模型一：

进行正向仿真设计光线追迹[4]

主要思路：

1．计算反射光强度

我们把线光源和抛物面都离散化成为离散点，针对离散后的线光源上的每一个点，与抛物面上的各个离散点，可以求出入射光线。

由于我们已知了抛物面的曲面方程，也就可以求得该抛物面上点的法向量。

根据反射定律，由下面的公式

（1），根据入射光线和法向量可以求得反射光线。

公式

（1）『1』

推导和证明见附录1

然后，由反射光线直线和25米处的测试屏相交，可求得交点坐标。

用双重循环，遍历线光源上和抛物面上的所有点，记录测试屏上通过B,C点（在误差应许的范围内，也可以是Ｂ，Ｃ附近的点）处光线的数量N。

２．计算直射光强度

为使直射光线与反射光线含同样的光强量。

（也就是跟反射的光线统一单位）。

只有这样才能考虑叠加效果。

只要把直射的光强也按反射光线同样的离散度转化为光线，这样就能统一单位了。

而反射光线是根据抛物面的离散程度来转化的，所以直射光线也需要经过抛物面来转化。

转化方法如下：

把测试屏上的B（或C）点与抛物面上的各个离散点直接相连。

得到直线方程，与线光源线段去交，（方法同上）

若有交点,则认为线光源上有一份光线对B（或C）点的光强有贡献。

若无,检测抛物面上的下一点,重复上述步骤.直到遍历完抛物线上所有点。

（记录个数M（直射））

３．求出最终结果

所得到的经过B点（或C）的光线数的和M+N，就可以表示线光源直射和反射到B（或C）点的光强。

同样,求出C点光强,叠加后,比较B与C点的光强之比。

重复上述步骤１，２，３，线光源线段（15,y,0）范围线光源按步长1——>0.001由大变小，逐步搜索，最终得出最优解为3.9280。

计算机实现[2]：

1．将抛物反射面离散化，划分为m个面积相同的网格，取得每个网格中心点的空间坐标数据。

2．线光源从一定步长d增长，没增加一个步长，将此时长度的线光源离散化为n个点。

取的它的空间坐标数据。

3．对线光源上的离散点和抛物面上的离散点进行两两连接，这样就取得此时长度下线光源发出的射向抛物面的所有入射光线。

4．利用MATLAB上的函数surfnorm求得抛物面（曲面）上每个离散点的法线向量。

5．根据入射光线和法线向量，也就是相当于已知入射光线向量和法线向量，通过公式

（1），我们可以求得每条入射光线通过抛物面反射面反射后的反射光线。

6．根据已知条件，测试屏上的两点B（25000，1300，0）和C（25000，2600，0），所以我们在空间上取x＝25000这个平面做为测试屏。

对5步鄹所求得的所有反射光线求取其在x＝25000平面上的交点，记录下所有点的坐标数据。

7．以B点和C点为圆心，在平面x＝25000上做一微小半径r的圆，记录在此小圆内的通过6步鄹求得的交点的个数，分别做为此时长度线光源通过反射在B点和C点产生的光强度（即反射光强度）的衡量值Bf和Cf。

8．将抛物面上的所有离散点和B点以及C点连线，分别计算连接B点的直线与线光源线段的交点个数和连接C点的直线与线光源线段的交点个数，分别做为此时长度线光源通过直射在B点和C点产生的光强度（直射光强度）的衡量值Bz和Cz。

9．因为线光源长度是以一定的步长增加的，所以可以取得满足条件（Bf＋Bz）>=2（Cf+Cz）（即B点光强不小于C点光强的两倍）的一系列长度值。

取得长度最大的一个长度值做为线光源的最优长度值，此时所需的功率是满足条件（Bf＋Bz）>=2（Cf+Cz）下的最小功率。

10．解过程中的数据图标如下：

B点光强

627

C点光强

292

297

强度比

2.14726

2.11111

线光源长度

3.91

3.912

3.914

3.916

3.918

627

632

297

302

2.11111

2.12795

2.09272

3.92

3.922

3.924

3.926

3.928

11．根据以上数据和算法原理，我们得出线光源长度的最优解为3.9280mm。

此时B点的光强度和C点的光强度的比值为2.09272。

（源程序见附录2）

模型二：

逆向光线追迹法

主要思路：

根据光路可逆原理，从B，C点出发进行逆向光线追迹，即由反射光线去寻找入射光线。

入射光线与线光源之间的公共点的存在，即表示所考虑的光线对B，C点的光强有贡献。

再由B，C点光强数（光线数）和倍数关系，以及功率越小越好的约束条件，线光源逐步长搜索最优长度。

详细步骤：

对于某一定长的线光源，因为线光源上的各个的点光源的光照强度一样，那么，在点Ｂ，Ｃ处的光照强度的比值就等于分别对于Ｂ，Ｃ点光源有贡献的光线数的比值，分析光线射到B点的照射情况：

１．考虑反射效果

由B（或C）点与抛物面上的离散点Ｍ（x0,y0,z0　），那么就可以求出入射光线了。

由于抛物面的方程：

F=z2＋y2-60x，

那么，在M（x0,y0,z0）处的法向量为：

（FX，FY，FZ）=（-30，y0,z0）

由计算原理中的公式（１）

可求出反射光线,该反射光线方程与线光源所在的直线去交。

并且，y取值范围就是线光源的长度.若有交点,则认为光线对B（或C）处的光强有贡献。

遍历抛物面上的所有离散点，记录下这些对Ｂ（或Ｃ）有贡献的光线的个数。

（记录个数N（反射））

２．考虑直射效果

为使直射光线与反射光线含同样的光强量。

（也就是跟反射的光线统一单位）。

只有这样才能考虑叠加效果。

只要把直射的光强也按反射光线同样的离散度转化为光线，这样就能统一单位了。

而反射光线是根据抛物面的离散程度来转化的，所以直射光线也需要经过抛物面来转化。

转化方法如下：

把测试屏上的B（或C）点与抛物面上的各个离散点直接相连。

得到直线方程，与线光源线段去交，（方法同上）

若有交点,则认为线光源上有一份光线对B（或C）点的光强有贡献。

若无,检测抛物面上的下一点,重复上述步骤.直到遍历完抛物线上所有点。

（记录个数M（直射））

３．求出最终结果

所得到的经过B点（或C）的光线数的和M+N，就可以表示线光源直射和反射到B（或C）点的光强。

同样,求出C点光强,叠加后,比较B与C点的光强之比。

然后，重复上述步骤１，２，３线光源线段（15,y,0）范围线光源长由大变小，逐步搜索，比较.最后可得出最优长度为４.06（mm）。

模型二的计算机实现[2]：

1．将反射面（抛物面）分解为721

721个等面积的网格，得到721

721个网格中心点的坐标数据集合X[721,721]，Y[721,721]，Z[721,721]。

2．利用MATLAB中的函数surfnorm计算得到反射面（抛物面）上每个网格中心点对应的法线单位向量NX[721,721],NY[721,721],NZ[721,721]。

3．在空间中取定测试屏上的两个测试点B点和C点（B（25000,1300,0）,C（25000,2600,0））。

4．分别从B点和C点向反射面（抛物面）上的每个网格中心点引向量。

经过B点的向量数组为（Bx[721,721],By[721,721],Bz[721,721]），经过C点的向量数组为（Cx[721,721],Cy[721,721],Cz[721,721]）。

即从线光源发出的光线镜抛物反射面反射后经过B点与C点的反射光线向量B，C.这样就可以利用每个向量和对应网格中心点的空间坐标数据唯一的确定每条反射光线。

5．利用公式求出每条经过B点和C点的反射光线的入射光线向量（Brx[721,721],Bry[721,721],Brz[721,721]）和（Crx[721,721],Cry[721,721],Crz[721,721]）

这样我们同样可以利用每个入射光线向量和对应网格中心电表数据确定每一条入射光线。

6．分别对经过反射面反射到B点的入射光线和经过反射面反射到C点的入射光线求与线光源所在直线的交点个数，分别记为Mfb和Mfc，作为线光源在此长度下对B点和C点的反射光强度衡量值。

7．分别求出经过B点和C点的反射光线与线光源所在直线的交点个数，并分别记为线光源在此长度下对B点和C点的直射光强度衡量值Mzb和Mzc。

8．线光源长度以1——>0.001的步长从0开始增长，用按步长搜索法求出长度最长的同时又满足（Mfb+Mzb）>=2（Mfc+Mzc）的线光源长度值（我们认为长度越长功率越小）。

9．求解过程中的数据表格如下：

半长度

1.9

1.91

1.92

1.93

1.94

1.95

1.96

1.97

b点光线数

309

310

311

313

c点光线数

134

137

139

140

143

144

146

光强比

2.30597

2.26277

2.23022

2.21429

2.17483

2.15972

2.14384

半长度

1.98

1.99

2.01

2.02

2.03

2.04

2.05

b点光线数

313

317

c点光线数

146

149

154

155

158

159

161

光强比

2.14384

2.10067

2.05844

2.04516

2.00633

1.99371

1.96894

半长度

2.06

2.07

2.08

2.09

2.1

2.11

2.12

2.13

b点光线数

317

319

321

323

c点光线数

161

162

164

166

169

173

174

光强比

1.96894

1.95679

1.93293

1.92169

1.89941

1.86705

1.85632

半长度

2.14

2.15

2.16

2.17

2.18

2.19

2.2

b点光线数

323

325

c点光线数

176

179

180

181

光强比

1.83523

1.81564

1.80556

1.79558

10．根据以上数据和算法原理，我们得出线光源长度的最优解为4.0600mm。

此时B点的光强度和C点的光强度的比值为2.00633。

（源程序见附录3）

模型三　:

具体思路

１．考虑反射的情况：

反向光束法是用来求由光反射引起测试屏上的点B,C处的照度大小的情况（光线追踪与模型二的情况一样）.反向光束法的思想包括沿着真实光线相反的方向追踪一个想象的光束。

射向照明系统的反向光束与光学系统的部件相交（抛物面反射镜）。

这些交点可以等效地看成发光源，称为等价光源。

它假想发射的光线到点B,C处照度的影响与点光源经过反射到点B,C处产生的照度影响等价（照度虽然与光源到点B,C的距离的平方成反比,但是由于车灯到测试屏的距离很大,在车灯里的点光源由于位置的不同所引起的与点B,C的距离不同的因素就忽略不计）。

２．考虑直射的情况。

对于某一个线光源上的点光源,它直射引起的点光源个数为1（即线光源上的该点光源）,３。

由于Ｂ，Ｃ照度是有同一个线光源产生的。

又由于线光源上的每个点光源对Ｂ，Ｃ的照度的影响是一样的。

则，点Ｂ，Ｃ处的照度也就与点光源（包括等价光源）的个数成正比.那么,点Ｂ，Ｃ的照度比就是可以反射光线到Ｂ，Ｃ的点光源（包括等价光源）的个数比。

.具体的做法如下:

先表示线光源在点B处时的照度情况

1）:

在抛物面上设有一点M（x0,y0,z0）,那么到点B（25015,1300,0）（单位:

毫米）的反射光线的向量为（x0-25015,y0-1300,z0）.由于抛物面的方程为：

F=z2＋y2-60x

那么在M（x0,y0,z0）处的法向量为：

（FX，FY，FZ）=（-30，y0,z0）

设入射向量

为（Rx，Ry，Rz）那么由公式（１）

求得：

Rx=-1799x0+1799*25015+60y02-60*1300y0+60z2

（２）

Ry=60x0y0–60*25015y0+y0-1300-2y03+2600y02-2y0z02（３）

Rz=60x0z0–60*25015z0+60y02-60*1300y0-2z03+z0（４）

2）:

这条入射光线的方程如下：

由于这条入射光线是由线光源上的点光源产生的，所以这条入射光线与直线x=15

z=0相交。

把x=15，z=0代入上面的式子，可以得出如下方程

（５）

3）又因为点M（x0,y0,z0）是在抛物面上，所以

y02+z02=60x0（６）

以上由式子（２）——（６）就可以得出y0　＝f（x0）,z0=g（x0）

4）有

，那么，与线光源z=0,x=15的交点处的y值可以用下表示

（0-z0）+y0=q（x0）（７）

从焦点开始，以一定的步长增加，对于每个y值yi，根据上面的六个方程，六个变量，就可以算出符合条件的所有点M（x0,y0,z0），这个点的个数就是等价光源的点的个数。

再加上一个直射的光源点个数1，设总个数为mi.

同理，

可以求出反射光线经过点C时候的情况，对应与每个与到点Ｂ情况时同样的y值yi，都有总的光源数（包括等价点光源）ni,那么根据题目的要求，使得线光源在点B，C处的光照强度符合1：

2的关系，只要逐步地以一定步长增加y的值，只要符合如下条件（取了y的个数为N）：

（

）/（250152+13002）

2（

）/（250152+26002）

记录这个时候的所有符合条件的yN。

在这些值里找yN　最大的情况

那么，２yN　也就是线光源的长度。

此模型是前面两个模型的改进，理论上会有更加精确的解，由于时间有限，我们未来得及用程序实现。

五、模型的评价和改进方法：

模型一的评价：

这种方法思路清晰,模型简单且利用计算机编程容易实现,误差取决于抛物面和线光源的离散精度.步长越小,精度越高,误差越小,误差分析简单.

但这种模型计算量随着精度的提高会成倍增加.在程序实现过程中,抛物面离散点数为m,线光源离散点数为n,则用循环计算入射光线、出射光线的时间复杂度就为0（m*n）。

而要得到较高的精度，我们一般需要把m设为50000以上，n设为100以上，那么时间复杂度至少为5.0e+6，因此此中分析算法的缺点是计算机处理时速度慢，效率低，需要较长的时间来处理。

此模型不仅适用于反射面为抛物面,还适用于椭圆面,反射镜,复合反射面反射镜,甚至曲面反射镜）。

对光源方面而言，此模型不仅适用于光源是线光源的光源设计，也适用于多线光源自由摆放的光源设计。

所以此模型适用性较广，但在高精度要求下，计算机运行速度必须相当快，才可能较快的得出理想的数据。

所以此模型在推广过程中较依赖于运行计算机的硬件环境。

模型二的评价：

逆向光线追迹法只对抛物面进行了点的离散，然后直接让利用通过B、C点的反射光线求得的入射光线与线光源直线段去相交。

这样做既避免了线光源的离散，提高了计算的精确度，又提高了计算机实现的程序的执行效率，因为循环次数的成倍减少大大缩短了程序的执行时间。

另外，逆向光线追迹法只针对测试屏上的B点和C点进行考虑。

题目的要求：

衡量线光源长度的取值是否合理，满足条件之一就是B、C点上光强度的比值关系，这使我们这样的做法成为合理和可行。

这样做，我们就避免了对一些无用点的计算，大大地减少了所需统计和计算的数据，而又不降低反而是提高了精确度。

我们尝试让用正向光线追迹的算法程序

（1）和逆向光线追迹的算法

（2）程序同时执行求解满足题意要求的线光源长度的任务，程序

（1）执行了将近1个半小时，而程序

（2）只执行了不到2分钟的时间。

模型三的评价：

这个模型中，使得点光源能够通过反射到点Ｂ，Ｃ处时，对于抛物面上的反射点M（x0,y0,z0）是有条件限制的。

这个限制可由式子（２）——（６）得出y0　＝f（x0）,

z0=g（x0），而模型二里只是离散地找了一些点，就会漏掉其它的同样符合条件的反射点。

所以，模型三是比模型二更加地精确，误差越小。

六、在有标尺的坐标系中画出测试屏上反射光的亮区的实现：

1、MATLAB软件实现：

取10m×10m的测试屏，将测试屏划分为41×41个方格，取得每个方格中心点的坐标数据。

对每个方格中心点利用逆向光线追迹的方法求得每个方格中心点的反射光强度（衡量值为光线条数），将光强度（光线条数）存贮为一个二维矩阵，这个二维矩阵就相当于一幅图象的颜色映射表，利用MATLAB里面的pcolor函数绘出该颜色映射表所对应的图象，此图象就是测试屏上反射光强度的分布图，颜色越浅表明光强度越大。

（源程序见附录4）

2、BorlandDelphi实现

在确定了线光源的长度以后，我们以如下的方法来寻找反射光的亮区:

1、在线光源上取一段微小的长度，并将其看作点光源S，其坐标为（Sx，Sy，Sz）；

2、把反射抛物面看作数量极大的反射点组成的集合；

3、按一定的顺序遍历所有的反射点，在遍历的过程中执行以下操作：

A设每次得到的反射点坐标为（X0,Y0,Z0），由抛物面方程可以得到该点的切平面（也就是相当于反射时的镜面）的法向量为（X0/30，Y0/30，-1）

B求出光源S关于切平面的对称点,也就是求出S的虚像。

这可以通过先求出直线L的方程，再得到L与切平面焦点A，通过点的对称来得到S’点的坐标。

C从S’点和反射点（X0,Y0,Z0），就可以确定反射光线所在直线M方程；

D光屏所在平面方程为Z＝25015mm,将其代入直线M的方程就可以得到反射光与屏的交点。

E大量交点的合成就是所求的亮区。

按上述的方法找到线光源上所有点对于反射抛物面上的所有点的反射光线与屏的交点。

程序运行结果为下图。

其中颜色越亮，表明该点的亮度值越大。

展开阅读全文