题车灯线光源的最优设计数学建模

1、A题 车灯线光源的最优设计参赛队员：王之元 谷德峰 饶彬 指导老师：毛紫阳学校：湖南长沙国防科技大学A题 车灯线光源的最优设计摘 要 车灯线光源的设计具有很强的实际应用意义。该问题属于单目标规划中的非线性规划问题。本文通过已知条件求出了灯光焦点，以及任一条反射光线的空间解析表达式和对应屏上的坐标位置表达式。然后建立光子跟踪模型进行求解。光子跟踪模型的原理是把光线粒子化，及时跟踪光子的运行方向，最后以单位面积打到屏上的光子数来衡量光照度大小，进而反映光强度在屏上的分布规律。这是一种离散型处理方法，其本质是计算机模拟。这个模型中基于不同原理又提出了好几种算法：等间距光子跟踪算法，改进的等间距跟踪算

2、法，等效立体角跟踪算法和随机方向跟踪算法。每一种算法的原理都不一样，层层递增，一步比一步深入，并分别作图进行比较。另外还结合边界条件讨论了线光源长度的临界值以及B、C两点光强随线光源长度变化的规律。在进一步讨论中，我们分析了光线直射到屏上时的情况，对二次反射的影响也做了分析，进而证明了问题的合理性。 本文我们得出的结论是：满足功率最小时的灯丝长度为4.337mm，第二问的答案见下图，图的大概形状是一个心形。关键词：光子跟踪模型 计算机模拟 一、问题重述与分析：1、 问题重述车灯线光源的设计是一个非常实际的问题。已知车灯的形状为一旋转抛物面，其开口半径36毫米，深度21.6毫米。经过车灯的焦点，

3、在与对称轴相垂直的水平方向，对称地放置一定长度的均匀分布的线光源。要求在某一设计规范标准下确定使光源功率最小的线光源长度（规范化要求略，见原题）；并对得到的线光源长度，在有标尺的坐标系中画出测试屏上反射光的亮区和讨论该设计规范的合理性。2、问题分析显然在线光源单位长度光通量一定的情况时，要使光源功率最小，线光源的长度也应该较小。但线光源的长度太小了，有可能出现C点的光强度小于额定值；线光源的长度过大，虽然能同时满足B、C两点光强度的要求，但线光源的功率也增大了。我们的目的就是在B、C两点光强度满足题目要求的情况下，求出最小的线光源长度。另外还要特别注意对“光强度”这一概念的理解，我们认为它和物

4、理学上的“发光强度”是一致的。按光度学中的定义，发光强度是某一方向上单位立体角内所辐射的光通量大小。一般不是用肉眼可以观察到的，主要的测量仪器是前照灯检测仪（参看中国汽车检测网前照灯检测）。其构件一般是采用具有把吸收的光能变成电流的光电池元件，按照前照灯主光轴照射光电池产生电流的比例，来测量前照灯的发光强度。由于本题中的光源不是点光源，直接求光强度比较困难，我们通过对光照度的测量来近似反映B、C两点光强度的大小。光照度是单位面积上所接受的光通量的大小，用来衡量被照明表面明暗程度的物理量。二、 模型假设：1、基本假设（1） 光线通过车灯的前玻璃时能量无损失。（2） 光线在抛物面进行一次反射时，能

5、量考虑成无损失。（3） 不考虑二次反射。（4） 不考虑光的干涉和衍射现象。（5） 截取线光源上很小的一段，可以看成是在空间呈均匀辐射的点光源。于是线光源可以看成是无穷多个点光源的叠加。（6） 不考虑线光源对反射光能量的阻挡和吸收，即是说线光源不考虑厚度，反射光线可以毫无影响地穿过线光源区域。2、 符号说明-车灯的深度，h=21.6mm-车灯的焦距-线光源的长度-单位长度上线光源的光通量，一个点光源的光通量可以近似记为- 一个点光源的光通量大小=- 一个点光源的发光强度 ，三、 模型建立及求解：定理1：灯泡的焦距=15mm证明：我们以光轴的正方向为x轴，竖直 方向为z轴，水平方向(AC方向)为y

6、轴，抛物面的顶点为原点建立空间直角坐标系。如图是抛物面在z=0时的函数图象，设函 数为则焦距= ，由抛物线的定义，得 解得p=30,所以=15mm。 #因此曲面方程：， 作出的图象见右图。该模型是一个单目标非线性规划问题其中K是某一要求的额定光强度, ,分别表示线光源在B,C两点的光强度 我们的目的就是求出满足P最小的 光子跟踪模型：本题中线光源可以看成是许多个点光源的叠加，只有焦点处的那个点光源通过抛物面后的反射光线是平行光束。其它点光源的反射光线都不是平行光束，因此射到屏上的总光线分布比较复杂，难以求出解析表达式。但并不是无规律所寻，给定了空间的一个点光源G(15,y1,0)，（其中 ），

7、再给定抛物面T: （上的一任点P(x0,y0,z0),则反射光线就可以根据y1,x0,y0,z0这四个参数唯一确定了。下面我们就根据光的微观粒子性结合反射光线的曲线分布建立光子跟踪模型。这一模型的主要思想是把光能粒子化，及时纪录每一个“光子”的运动轨迹，判断最终到达屏上的具体位置。最后把每个小面积区域内的“光子”总数进行叠加，就得到了屏上光通量的分布图，进而求出光照度和光强度的分布。例如对于点光源G(15,y1,0)，（ ），设其光通量为，单位时间内辐射出的“光子”数为M(为提高精度，M可以取得大一些，如1000，000个)，则每个光子所带能量为/M（一秒钟内）,且这些“光子”围绕G点在空间是

8、均匀发射的。一部分“光子”通过灯泡外玻璃壳直接射出（这一部分能量非常小，我们在模型的进一步讨论中将继续讨论），绝大多数的“光子”通过抛物面反射一次后射出，还有可能出现“光子”在抛物面反射两次，在反射第二次时，我们设想该“光子”就被牢牢的地粘在了抛物面上（这种情况概率也很小，模型改进中将继续讨论）。通过把线光源选取一定步长dl做为一个点光源，可以得到N个点光源(N=),我们把这N个点光源在屏上某一小区域内的光子数进行求和,就可以求得该小区域内的光照度值我们先分析一下任一个点光源G(15,y1,0)，（）在抛物面上任一点P(x0,y0,z0)的反射曲线.为了直观的理解我们在平面上做出了大致反射光线

9、图,实际上这些点图并不是在一个平面上的（见图2）。G是点光源，P是反射点，E是G点关于过P点切平面的镜面对称点，PN是反射方向。由空间解析几何的知识，我们不难得出E点坐标：进而由反射光线PN和屏面方程x=25015可以求出反射光线打在屏上的坐标K(x3,y3,z3),由Mathematica计算结果见下表。公式较长，但我们在用计算机进行模拟时需要用到数据采样：下面我们就给出具体做法，首先要进行数据采样，即是说给定一个点光源G，让它遍历抛物面上的所有点。我们遍历的参数是和 (,0,360degree)。 其中是抛物面上的一点到x轴的垂直距离，垂心为是以为圆心，为半径上的圆转过的角度。则抛物面方程

10、化为 （如图3）当我们的，定了时，抛物面上一点也定了。，分别变动很小的时，光线将在抛物面上扫出一个小面积dS1，这个面积非常小，我们可以近似的认为是一个小平面，光线照射就可以看成是镜面反射。我们设想有一个虚拟光源G, G 是G关于dS1的镜面对称点.则从G点反射出的光线完全可以等效看成是从G点直射出的光。设G点通过dS1射到屏上的面积为dS2,由光度学的知识可以得到dS2上的光照度其中dF是点光源G在dS1面上的光通量,d是G到dS2的距离，是从G射出光线与x轴的夹角。给定点光源G，给定，和 时，光线在屏上的位置和照度就唯一确定了。下面我们就建立算法来实现遍历过程。1、 等间距光子跟踪算法st

11、ep 1 选定点光源步长dl,确定点光源个数N=,每个点光源1秒内通过的光子数为设为M=100000,step 2 将吸收屏均匀分割成面积为ds许多小正方形,记下位置和编号step 3 选定,分别确定和的遍历次数, ,每次遍历对应的dS1通过的光子数记为相同step 4 for(I=1;I<=N;I+)for(j=1;j<=N1;j+)for(k=1;k<=N2;k+) 计算打到屏上的光子位置,累加每个小正方形上光子数, step 5 是否完毕,返回根据上述算法我们利用所得离散数据作出了光线在屏上分布的形状图, 上图是线光源长度等于6mm时做出的图,由图可以看出反射光线在屏上

12、的形状是心型。其它长度作出图的形状也大致呈心型，只有处在焦点处的那个点光源在屏上的光照度是大致呈圆形分布的（即是说线光源非常短，浓缩到一点的极端情况，0mm）（见图6）。这和实际情况也是比较符合的，而且可以看出光照度的分布规律是中间大（非常大），周围小 。但算法1有严重的理论不足之处，它只能反映反射光线在屏上的形状，而不能很好的反映在屏上的照度变化情况。出现这一情况的原因是在进行步长遍历时，把每一次对应的dS1上通过的光子数考虑是一样的了，实际上是不相同的。另外由于光屏离的很远（相对于抛物面的大小来说），在抛物面上的dS1很小，但发散到屏上时的光斑dS2可能比较大，如果还把它看成一个质点来处理

13、也是不合理的。光子从点光源向空间各个方向是均匀发射的，我们应该以通过等光子数的小单元进行遍历。这在以点光源为中心的球坐标上可以等角度遍历，但反映到以进行遍历时，通过相同光子数对应应该是变步长的，为了便于操作，我们仍然实行等步长遍历，只是对于每一步的光照度我们定义了一个权，这个权值的大小和，是有关系的，而且由于第二个原因的影响，故我们采用分布叠加的方式求屏上的光照度。（见图7）它的主要原理是一条光线反射到抛物面上后计算其光照度，然后又把这个反射点近似的看成是点光源，只不过是光强度变小了，再根据光照度公式可以求出屏上的光照度分布。这种算法很好的减少了由于距离太长带来的模型误差。 设点光源G(15,

14、y1,0)发出的一条光线与抛物面交于P,则其光照度,从p点反射后,我们可以把p点看成是一个光源,则打到屏上M点的光照度就为,再将每个点光源经抛物面反射的光线在屏上的光照度进行叠加就可以得出B、C点光照度的比例关系。有关计算结果如下，入射方向n1 : 出射方向i： Q即是我们要设定的权值，它和 有关，则照到屏上M点的光照度为然后我们对，等步长进行遍历，每一补加上相应的权值就可以得到光照度在屏上的分布规律。2、改进的等间距跟踪算法只需要把算法1的step 3改为step 3 计算等步长内反射到屏上加权光 照度由算法2我们得到的光照度在屏上的一组分布分布图（图8A），图8A是线光源长度=3.0mm时

15、的照度分布图。该图象从直观上看能很好的反映光的明暗情况。但这里得到的图象并不是我们在实际生活中所看到的光斑。光度学中唯一能用眼睛感受的物理量是光亮度，即是我们平常所说的亮度。但光照度的分布和光亮度的分布大致是相同的，可以用照度来反映亮度的分布情况。 3 等效立体角跟踪算法算法2虽然考虑到加权的因素，但加权是利用光强加权，显的不是很合理。对此我们提出等效立体角跟踪算法。其原理是变步长进行的遍历，但每一步所对应小面积上的立体角是相等的。每次遍历到下一点后的值和， 以及前一步的是有关的。我们找出了它们之间的关系就可以看成是等效相同的立体角进行遍历。图9是这个算法的流程图，图9A是这个算法得到的一组屏

16、上光照度分布图。,图 94 随机方向跟踪算法前面三种算法都是确定性算法，需要在抛物面上取样跟踪，而且还要等效处理变步长问题，计算比较复杂，而且容易产生模型误差。为此我们提出随机方向跟踪算法，这是一种概率算法。原理很简单，我们定义一个随机方向Pa,b,c,(a,b,c-1,1 ) 。利用计算机按均匀分布产生M对a,b,c的值，M一般都取的比较大（例如M=1000000）。只要M足够大，随机方向P就可以看成是在空间呈均匀散开的。我们只需要对每一个方向求出在屏上的照射点。累积屏上单位面积的反射光线条数，就可以得到光照度的分布。该方案执行简单，但也有不足之处，实际上空间分布应该呈球状均匀分布，而不是呈

17、立方体。我们也可以等效加权处理。步骤如下（1） 对于产生的随机方向P,如果落在以1为半径的球内,则权值为1.（2） 对于产生的随机方向P,如果落在以1为半径的球外,则权值为1/d.其中d为该点到球心的距离。运用算法3对程序进行少量改动，我们作出了屏上光照度的分布图（图8C）。和图9A相比两者几乎完全相同，说明我们的随机方向算法算法是可靠的。上面的处理方法也是一种近似方法，可以用数学论证严格的推出权值应满足的规律。另外由于是随机产生的发射方向，总有可能出现方向分布不是很均匀的情况，最好是进行方差检验或是增加划分细度。灯丝长度的范围的讨论：我们先从理论上来分析灯丝长度应满足的关系式。线光源发出的光

18、经抛物面反射后如果能直接照到B点和C点，则反射点应该在z=0的抛物线上。因此我们可以只研究z=0的抛物线上的反射规律来确定线光源长度的范围。线光源长度如果太短了，反射出的光可能照不到C点（极端情况线光源长度为0，全部反射出平行光，不能照到B,C点，随着线光源长度的增加，能逐渐照到B,C点。当线光源长度继续增加到某一值后，超过这一值的长度部分发出的光经反射后将不能照到C点，而会出现在C点以外。这一部分光线做的就是额外功，因此线光源的最优长度应在某两个特定植之间。我们在Mathematia中求得了给定点光源G(15,y1,0)情况下，反射光线斜率r随x0的变化函数其中x0是抛物线上的点，x0,21

19、.6,当y1给定了时r随x0的函数图象是单调变化的，图9给出的是y1=3mm的情况（反射点纵坐标y0>0的图象，）。由图可以看出斜率是单调下降的，x0=0处斜率最小，但不为0 ，原因是y1不等于0，y1越大，图象在y轴的截距也将大些。由图象的单调性我们可以推知，线光源长度下界应满足的临界条件是点光源G(15,y1,0)在（21.6,-36,0）处的反射光线恰好过C点，由于只考虑了（21.6,-36,0）的情况，还应该考虑（21.6,36,0）点的情况，二者求出的较小一个2 y1值就是长度下界。对于长度上界满足的邻界条件刚好相反。临界条件是点光源G(15,y1,0)在（21.6,36,0）

20、处的反射光线恰好过C点。经过解方程计算得到=/2 <=1.58202, 即3.07762<<3.16404 这就是灯丝长度的范围。 上面考虑的是单位长度光通量一定时候，线光源长度应满足的条件。如果我们把线光源的总功率看成是一定的，则单位长度的光通量是随线光源的长度而变化。因此的上界就不能用上面的方法去判定。但下界是一样的，因为下界是一个突变的过程，如果线光源长度少于下界值 =3.07762，是不可能有光到达C点的。 对于上界也可以定性判断规律，当线光源长度超过3.16404时，B,C两点的光强度是随的增大而减少的。这是因为超过3.16404的部分不会再照到B,C上，但由于总功

21、率一定，长度在3.07762<<3.16404中的点光源的光通量却减少了，因此引起光强度的下降。线光源最优长度值的确定：我们利用算法2得到的数据分析了线光源长度 的变化对B、C两点光强度影响的情况，拟合曲线见图10。 图中x轴表式灯丝长度，y轴表示光照度，由于我们采用的是光子跟踪模型，y轴的值比较大，但并不影响B,C两点光照度的相对关系。最优线光源长度只和他们之间光照度之比有关。如果把线光源的功率给定了（光通量也定了），则线光源的长度越长，单位长度的光通量越小，到达B,C两点的光强度越小。极端情况线光源无穷长，超过临界长度的发光部分都浪费了，射到B,C点的光强度就十分小，不能满足题

22、目要求的达到某一额定值。直接从该图上还不能确定最优线光源的长度值。所以我们作出了在线光源功率一定的情况下B,C点光照度之比随线光源长度变化的函数图象 。 当B点光照度小于等于C点光照度2倍时,制约满足原题条件的额定功率P大小的必要条件是B点光照度,因为,当B点光照度满足>=2P时,C点光照度必满足>=P;当B点光照度大于等于C点光照度2倍时, 制约额定功率P大小的必要条件是C点光照度,因为,当C点光照度满足>=P时,B点光照度必满足>=2P;我们是在功率一定的条件下对线光源长度进行模拟搜索,所得数据点经过拟合处理后,得到上面的两条曲线,当L增大时;起初, B点光照度大于

23、等于C点光照度2倍, B是制约额定功率P大小的主要因素,经过一段距离后,B点光照度小于等于C点光照度2倍, C成为制约额定功率P大小的主要因素;综合分析得L/2=2.1685倍附近时,得到满足在功率一定的条件下的最大额定功率条件,而原题恰好是一个逆向思维,当额定功率一定时,分析求解满足最小功率条件的线光源的长度,可见,这于固定功率,求可满足的最大额定功率是等价的。我们作出了在线光源长度L=4.337mm情况下屏上的反射光亮区（见图11），这即是第二问的答案四、 模型的进一步讨论： 前面的模型都是不考虑直射的情况得出来的，实际上射到屏上的一部分光还有来自于灯泡直射的部分。因此严格的说光斑的形状还

24、要大些，但也不是圆形，而应该接近于椭圆型。只是由于直射部分的能量占的比重太小，看起来还是心型起绝对作用。但这也反映了设计规范性的合理问题。 下面我们就来推导一下直射时光照度应满足的规律设点光源G(15,y1,0) ，所带光通量为Qdy1,直射到屏上AC方向的坐标为M(25015,y3,0)则该点光源在M点的光照度为 那么线光源在M点的光照度就是这些点光源在该点的积分(y1)将这一部分光强也考虑进去，我们得到更加合理的线光源长度值L=4.337mm这个值比仅考虑反射的情况要稍微大一些。但也可以看出直射光线的能量相对于反射光线的能量是很小的一部分，近似计算中完全可以不考虑。另外还需要考虑2次反射的情况，可以根据反射光线与抛物面方程的交点来