2002年数学建模A题——车灯线光源的优化设计

题目要求

注：问题三略

问题分析

问题一

在满足该设计规范的条件下，计算线光源长度，使线光源的功率最小。其中设计规范为：
要求C 点的光强度不小于某一额定值(可取为1 个单位)，B点的光强度不小于该额定值的两倍(只须考虑一次反射)。

显然，这是一个优化问题：

• 决策变量为线光源的长度 $l$ ；
• 目标函数为线光源的功率$W\to min$；
• 约束条件为：
  • C 点的光强度 $\ge$ 1(或K)
  • B 点的光强度 $\ge$ 2(或2K)

若用$W(l)$表示线光源的功率，$E_B(l)$ 和$E_C(l)$分别表示$B$和$C$点的光强度，则数学模型为：

   $minW(l)$
$s.t.$ $E_B(l)\ge 2K$
         ~~~~~~~~          $E_C(l)\ge K$

进而，需要解决的问题归结如下：

1. 如何表达目标函数$W(l)$
2. 如何表达约束条件，即如何表达测试屏上$B$和$C$点的光强度
3. 如何求解优化模型

问题二

当问题一求解出时，问题二可在问题一的基础上的更一般的问题。

建立模型

问题一

1. 目标函数$W(l)$

由于线光源的功率与长度正相关，因此当线光源长度最小时，功率便达到了最小。

2. 约束条件

通过文献可知，测试屏上某点的光强度应为光照度。
点光源在测试屏上某点的光照度为
$$E=\frac{Icos\theta }{R^2}$$
其中 $I$ 为点光源的发光强度，$\theta$ 为光源发出的光束其轴线与受照点处法线 $n$ 所成的角度，$R$ 为点光源到受照点
的距离。

根据题意，建立如下坐标系：

为了统一单位，我们把所有单位化为$mm$

1) 表示各点坐标

设抛物面方程为 $x^2+y^2=2pz$ ，将点 $G(36,0,21.6)$ 代入求得 $p$ 。

我们记线光源上的点为 $P(0,w,p/2)$ ，抛物面上反射点为 $Q(x,y,z)$ ，测试屏上点为 $T(x_0,y_0,25000+p/2)$
（设 $x^2+y^2=r^2$，则 $z=r^2/(2p)$）， MATLAB表示如下：

x=36; y=0; z=21.6;  p = (x^2+y^2)/2/z;
syms w x y r x0 y0 real,  z = r^2/(2*p);
P = [0, w, p/2];  Q = [x, y, z];  T = [x0, y0, 25000+p/2];

2) 表示各个向量

1. 求Q点处的法向量 n
   将抛物面记为$F(x,y,z)=x^2+y^2-2pz=0$，则$n=[F_x,F_y,F_z]$。
2. 还有向量PQ，向量QT。

这三个方向向量构成菱形及其对角线，如下图：

其中a为PQ，b为n，c为QT

syms z1,  F = x^2+y^2-2*p*z1;
F_x = diff(F,x);  F_y = diff(F,y);  F_z = diff(F,z1);  
n = [F_x, F_y, F_z];
PQ = Q-P;  QT = T-Q;

3) 寻找P、Q、T的关系

$$b=2\frac{a\cdot n}{|n|}\frac{n}{|n|},c=a-b$$

n = 2*dot(PQ,n)/dot(n,n)*n;    %//n = 2*dot(PQ,n)/sqrt(dot(n,n))*(n/sqrt(dot(n,n)))
QT1 = PQ-n;
syms k, clear F
F = QT1 - k*QT;      %//由于两向量平行，因此F=0
F = subs(F, x0, 0);  %//由于B点和C点x均为零，故令x0=0

经计算，$F_x,F_y,F_z$ 分别为：



人工化简得（注意$x^2+y^2=r^2$）：

$F_x=x\left(k+\frac{2wy}{r^2+90}\right)=0$

$F_y=ky+\frac{2w{y}^2}{r^2+900}-w-k{y}_0=0$

$F_z=\frac{r^2}{60}(k+1)-\frac{60wy}{r^2+900}-25015k+15=0$

• 当$x=0$时，$r^2=y^2$，$F_y,F_z$可进一步化简为：

  $F_y=k(y-y_0)+\frac{2w{y}^2}{y^2+900}-w=0$

  $F_z=\frac{y^2}{60}(k+1)-\frac{60wy}{y^2+900}-25015k+15=0$

  联立以上两式，消去$k$，化简，合并同类项：

  F(1) = 0;
  F(2) = k*(y-y0) + 2*w*y^2/(y^2+900) - w;
  F(3) = y^2/60*(k+1) - 60*w*y/(y^2+900) -25015*k + 15;
  k0 = solve(F(2),k);   %// k0=(w - (2*w*y^2)/(y^2 + 900))/(y - y0)
  F(3) = subs( F(3), k, (w - (2*w*y^2)/(y^2 + 900))/(y - y0) );
  [f1, ans] = numden(F(3));
  f1 = collect(f1, y);
  

  得到 $f_1=$

  

---

未完待续