【2021数模国赛 江苏省二等奖】基于“FAST“主动反射面的形状调节研究

基于"FAST"主动反射面的形状调节研究

摘要

本文利用几何原理，欧拉角，非线性规划与三向量共面以及计算机编程确定理想抛物面，进而通过变形策略确定工作抛物面以及调节后反射面300米口径内的主索节点编号、位置坐标、各促动器的伸缩量等，最后获得天体电磁波经反射面反射后的最佳接收效果。

针对问题一，当待观测天体 $S_1$ 位于基准球面正上方，即 $\alpha =0°$ ，$$\beta =90°$$ 时，并且由于基准球面与抛物线面都具有各向同性，可以将调节中的基准球面与抛物面的相关计算简化为二维平面下进行。变形策略主要考察径向位移最小，基准球面与抛物面方程转化为极坐标方程，结合考虑反射面板调节因素，再运用二重积分根据积分范围求各点位移之和最小值为0.0570以及最大伸缩量为0.3779，最终确定理想抛物面方程:

$$z=F-R-\frac{1}{4\times 0.001783589476658}+0.001783589476658(x^2+y^2)$$
( $R$ 为基准球面的半径， $F=0.466R$ )

针对问题二，借助欧拉角，将此时位于 $\alpha =36.795°$ ，$$\beta =78.169°$$ 观测天体 $S_2$ 以原点 $C$ 建立坐标系转换公式，利用问题一方法和理想抛物面方程，得出问题二的理想抛物面方程进一步数值计算得出问题二的理想抛物面没有超过基准球面的边界。在此基础上，先利用迭代法对主索节点进行调整，进而运用非线性规划进行优化，尽量贴近理想抛物面，并求得调节后反射面300米口径内的主索节点位置坐标等结果。

针对问题三，利空间几何向量法和物理光反射原理，求解出平行电磁波经反射面板反射，反射到馈源舱所在平面的点的坐标，得出面积。进而推广求解平行电磁波经反射面300米口径有效照明区域内反射面板总块数与圆盘重合累加总面积占总反射面积的比例，即馈源舱的接收比，等于0.060536495085120（约6.05%）。馈源舱的接受比小于基准反射球面的接收比0.069802419471494（约6.98%），是由于两者的每块反射面板接收平行电磁波面积不同。

关键词：FAST 二重积分 欧拉角 非线性规划 迭代法 三向量共面

一、问题的提出

中国天眼"FAST"，具有中国独立自主知识产权，是世界上目前口径最大、最具"威力"的单孔径射电望远镜，其设计综合体现了我国高技术创新能力。"FAST"作为多学科基础研究平台，有能力将我们的视野拓展至宇宙边缘，检测引力波，观测暗物质和暗能量，寻找地外文明和探究宇宙起源和演化。它在日地环境研究、国防建设和国家安全等方面发挥不可替代的作用。

它的500米口径射电望远镜反射面采用柔性索网作为主体支撑结构，通过促动器连接下拉索控制主索节点变位实现巡天观测。在进行天文观测的工作态时，对部分节点进一步拉紧而对另一些节点进行放松，从而在口径球面上形成一个焦径比为0.466，口径为300米的抛物面，在抛物面的焦点处放置一个馈源舱（接收装置），以对所接收的宇宙射电信号进行聚焦。通过控制系统使抛物面的焦点始终与馈源舱同步移动以便主动观测某天体。

在反射面板调节^[1]约束下，如何确定一个理想抛物面，然后通过调节促动器的径向伸缩量，如何将反射面调节为工作抛物面，使得该工作抛物面尽量贴近理想抛物面，以获得天体电磁波经反射面反射后的最佳接收效果成为迫切需要解决的问题。

需要建立模型讨论并解决下列的几个问题：

1. 当待观测天体 $S$ 位于基准球面正上方，即 $\alpha =0°$ ，$\beta =90°$ 时，结合考虑反射面板调节因素，确定理想抛物面。

2. 当待观测天体 $S$ 位于 $\alpha =36.795°$ ，$\beta =78.169°$ 时，确定理想抛物面。建立反射面板调节模型，调节相关促动器的伸缩量，使反射面尽量贴近该理想抛物面。将理想抛物面的顶点坐标，以及调节后反射面300米口径内的主索节点编号、位置坐标、各促动器的伸缩量等结果按照规定的格式（见附件4)保存在"result.xlsx"文件中。

3. 基于第2问的反射面调节方案，计算调节后馈源舱的接收比，即馈源舱有效区域接收到的反射信号与300米口径内反射面的反射信号之比，并与基准反射球面的接收比作比较。

图1.反射面板示意图

图2.主索网结构及其连接示意图

图3.FAST剖面示意图

二、基本假设

（1）索网结构不会因为风速，温度等环境影响对模型有任何影响

（2）不考虑下拉索与促动器的故障问题

三、符号说明

符号	说明
$S$	角度为 $\alpha$,$\beta$ 的观测天体
$R$	标准圆半径
$r$	标准圆口径
$r_1$	基准球面到$c$点的距离
$r_2$	理想抛物面到$c$点的距离
$r_3$	抛物线口径
$\theta$	有效照明角度
$R_j^i$	$j$系坐标转换到$i$系坐标的矩阵
$△M_1{M}_2{M}_3$	主索节点为$M_1,M_2,M_3$的反射面板
$S_{△M_1^{\prime }{M}_2^{\prime }{M}_3^{\prime }}$	平行电磁波经反射面板投影在馈源舱所在平面形成的面积
$k_1$	馈源舱的接收比
$k_2$	基准反射球面的接收比

四、问题1模型的分析、建立与求解

4.1 模型分析

根据当待观测天体 $S_1$ 位于基准球面正上方，即 $\alpha =0°$ ，$\beta =90°$ 时，由于基准球面与抛物线面都具有旋转对称性，调节形成的抛物面在球面上移动时与球面的相对关系保持不变。所以在分析反射面板调节时，可以将调节中的抛物面与球面的相关计算简化为二维平面下进行。以基准球面球心为坐标原点，建立基准球面和工作抛物面的方程。变形策略主要考察径向位移最小^[2]，即弧长$AB$上各点经过变形为$A^{\prime }{B}^{\prime }$所产生的位移最小值，将方程转化为极坐标方程，结合考虑反射面板调节因素，再运用二重积分求各点位移之和最小值以及最大伸缩量，最终确定理想抛物面。

4.2 模型建立

对问题一，由于基准球面与抛物线面都具有旋转对称性，故我们将抛物面与球面的相关计算简化到二维平面下进行。

由反射面为半径约300米、口径为500米的球面，可得基准球面公式：
$$x^2+y^2+z^2=300^2(-250\le x,y\le 250)$$

以z轴为对称轴的旋转抛物面的公式：
$$z=a+b(x^2+y^2)$$
其中，抛物面焦点为：$(0,0,a+\frac{1}{4b})$ ，P点坐标为：$(0,0,F-R)$ ，同时，由于P为抛物面焦点，所以有：$F-R=a+\frac{1}{4b}$ ，于是，可以得到旋转抛物面公式为：
$$z=F-R-\frac{1}{4b}+b(x^2+y^2)$$
此时，只有一个未知变量b，然后确定它的最优解，以此我们想到两种方式，一种是在径向变化最小的约束条件下进行，另一种是在弧长变化最小的约束条件下进行。题目交代是通过调节促动器的径向伸缩量，将反射面调节为工作抛物面，使得该工作抛物面尽量贴近理想抛物面，由此我们舍弃在弧长变化最小的约束条件下进行目标优化的方法，采用在径向变化最小的约束条件下进行。

图4.径向最小几何模型

由于标准球面与抛物线方程都具有旋转对称性，故考虑二维情况，二维平面模型如图4所示，即弧长$AB$上各点经过变形为$A^{\prime }{B}^{\prime }$所产生的位移最小值，将方程转化为极坐标方程，将位移化为二重积分求解。

4.3 模型求解

根据工作状态时，反射面的形状被调节为一个300米口径的近似旋转抛物面（工作抛物面），与反射面为半径$R=300m$的球面，再假设抛物面的顶端边界点和基准球面重合，易得全等三角形，则$\theta \in [-\frac{2}{3}\pi ,-\frac{1}{3}\pi ]$ ，此时直角坐标系下，基准球面方程为：$x^2+y^2=300^2(-250\le x\le 250)$ ，抛物线公式为：$y=F-R-\frac{1}{4b}+b{x}^2$ ；极坐标系下，基准球面方程为：$r1=300,\theta \in [-\pi ,\pi ]$ ，抛物线公式为：$b{r_2}^{2}co{s}^2\theta -r_{2}sin\theta +F-R-\frac{1}{4b}=0\theta \in [-\frac{2}{3}\pi ,-\frac{1}{3}\pi ]$ ，即：$r_2=\frac{\sin \theta +\sqrt{si{n}^2\theta -4b{\cos }^2\theta (F-R-\frac{1}{4b})}}{2b{\cos }^2\theta }$ ，化简可得：$r_z=\frac{300-\sin \theta +\sqrt{{\sin }^2\theta -4b{\cos }^2\theta (F-R)}}{2b{\cos }^2\theta }$ ,于是目标函数为：$\text{min}{\int }_{-\frac{\pi }{2\pi }}^{-\frac{\pi }{3}}(r_1-r_2{)}^{2}d\theta .$

即:
$$\begin{array}{c}min{\int }_{-\frac{2\pi }{3}}^{-\frac{\pi }{3}}(300-\frac{sin\theta +\sqrt{1-4b{\cos }^2\theta (F-R)}}{2b{\cos }^2\theta }{)}^{2}d\theta \\ \mathfrak{s}.\mathfrak{t}.|r_1-r_2\le 0.6|\end{array}$$

带入数据计算，得到b=0.001783589476658，于是，抛物面方程为：

$$z=F-R-\frac{1}{4\times 0.001783589476658}+0.001783589476658(x^2+y^2)$$

图5.理想抛物面方程

此时：

$$min{\int }_{-\frac{2\pi }{3}}^{-\frac{\pi }{3}}{\left(300-\frac{sin\theta +\sqrt{1-4b{\cos }^2\theta (F-R)}}{2b{\cos }^2\theta }\right)}^{2}d\theta =0.05700521145460$$

在$\theta \in [-\frac{2}{3}\pi ,-\frac{1}{3}\pi ]$上最大伸缩量为0.3779米。

图6.300m口径抛物面促动器行程内抛物面与球面偏离最大距离

五、问题2模型的分析、建立与求解

5.1 模型分析

与第一问不同的是此时观测天体$S_2$位于 $\alpha =36.795°$ ，$\beta =78.169°$ ，我们可以设$S_{1}C$为$Z_i$轴（以$\overrightarrow{C{Z}_i}$为正方向，建立右手坐标系），$S_{2}C$为$Z_j$轴（以$\overrightarrow{C{Z}_j}$为正方向建立右手坐标系），空间一点$P$从$j$系坐标利用欧拉角转换矩阵得到$i$系坐标，利用这一想法我们就可以清晰得出问题二所求理想抛物面的点都可以利用方向余弦转换到问题一的理想抛物面上，就可以得出__63.wmf)的理想抛物面，但还要分析问题二的理想抛物面是否超过基准球面的边界。

5.2 欧拉角转换

从一个坐标系到另一个坐标系的转换有多种方法：欧拉角^[2]法、方向余弦矩阵法、四元数法等。本文采取欧拉角，其核心思想是：一个坐标可以用另一个参考坐标系的三次空间旋转来表示。

设有两个共原点的右手系$O{X}_i{Y}_i{Z}_i$和$O{X}_j{Y}_j{Z}_j$，空间有一点P该点在两个坐标系内的坐标分别为$[x_i,y_i,z_i{]}^T$和$[x_j,y_j,z_j{]}^T$，P点从$j$系变换到$i$系的坐标变换为：

$$\begin{array}{c}\{\begin{array}{l}{x}_i=x_{j}cos(x_i,x_j)+y_{j}cos(x_i,x_j)+z_{j}cos(x_i,z_j)\\ y_i=x_{j}cos(y_i,x_j)+y_{j}cos(y_i,y_j)+z_{j}cos(y_i,z_j).\\ z_i=x_{j}cos(z_i,x_j)+y_{i}cos(z_i,y_j)+z_{j}cos(z_i{z}_j)\end{array}\\ [r{]}_i=R_j^i[r{]}_j\\ R_j^i=\left(\begin{array}{lll}cos(x_i,x_j)& cos(x_i,y_j)& cos(x_i,z_j)\\ cos(y_i,x_j)& cos(y_i,y_j)& cos(y_i,z_j)\\ cos(z_i,x_j)& cos(z_i,y_j)& cos(z_i,z_j)\end{array}\right)\end{array}$$

即为一般形式的转动矩阵，也称为从j系变换到i系的转动矩阵。对$i$系来说，$R_j^i$描述了$j$系的姿态，故也称其为姿态矩阵，又因$R_j^i$内各元素皆为坐标轴之间的方向余弦，所以也可称其为方向余弦矩阵。

当两个坐标系无相对转动时，$R_j^i=I$；下面求绕x轴的旋转的转动矩阵。

绕X，Y，Z的坐标轴的旋转变换矩阵是最基本的转动矩阵，他们是一般转动矩阵的特例，故可直接由一般转动矩阵得到。

图7.坐标系旋转示意图

可得到绕x轴旋转$\theta$角的转动矩阵为：

$$\begin{array}{c}R(X,\theta )=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\theta & -sin\theta \\ 0& sin\theta & cos\theta \end{array}\right]\\ R(Y,\theta )=\left[\begin{array}{lll}cos\theta & 0& sin\theta \\ 0& 1& 0\\ -sin\theta & 0& \cos \theta \end{array}\right]\\ R(Z,\theta )=\left[\begin{array}{lll}\cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]\end{array}$$

下面求绕两个坐标轴旋转的转动矩阵。

图8.内在旋转示意图

如图8所示，设坐标系$O{X}_i{Y}_i{Z}_i$先绕$Z_i$旋转$\theta$角形成坐标系$O{X}_{im}{Y}_{im}{Z}_{im}$，再绕$Y_m$轴旋转$\alpha$角，形成坐标系$O{X}_j{Y}_j{Z}_j$。

其中：
$$\begin{array}{rl}{\left[r\right]}_i& =[R_m^i(Z,\theta )][r{]}_m\\ {\left[r\right]}_m& =[R_j^m(Y,\alpha )][r{]}_j\\ {\left[r\right]}_i& =\left[R_m^i(Z,\theta )\right][r{]}_m=\left[R_m^i(Z,\theta )\right]\left[R_j^m(Y,\alpha )\right][r{]}_j\end{array}$$
依照定义，上述旋转的旋转轴是前面一次变换后的新系的轴，因此是内在旋转。可见，内在旋转要将绕各轴的旋转矩阵从左至右连乘。此时的旋转矩阵为：

$$\begin{array}{c}{R}_j^i(\theta ,\alpha )=\left[\begin{array}{lll}cos(\theta )& -sin(\theta )& 0\\ sin(\theta )& cos(\theta )& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{lll}cos(\alpha )& 0& sin(\alpha )\\ 0& 1& 0\\ -sin(\alpha )& 0& cos(\alpha )\end{array}\right]=\left[\begin{array}{lll}cos(\theta )cos(\alpha )& -sin(\theta )& cos(\theta )sin(\alpha )\\ sin(\theta )cos(\alpha )& cos(\theta )& sin(\theta )sin(\alpha )\\ -sin(\alpha )& 0& cos(\alpha )\end{array}\right]\end{array}$$

$$[r]=[R_m^i(Z,\theta )][r{]}_m\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right]=\left[\begin{array}{c}\cos \theta \cos \alpha x_j-\sin \theta y_j+\cos \theta \sin {\alpha }_j\\ \sin \theta \cos \alpha x_j+\cos \theta y_j+\sin \theta \sin \alpha z_j\\ -\sin \alpha x_j+\cos \alpha z_j\end{array}\right]$$

此时的旋转抛物面方程为：

$$\begin{array}{r}-\sin \alpha x_j+\cos \alpha z_j=F-R-\frac{1}{4b}+b\left[{\left(\cos \theta bs\alpha x_j-\sin \theta y_j+\cos \theta \sin \alpha x_j\right)}^2+{\left(\sin \theta \cos \alpha x_j+\cos \theta y_j+\sin \theta \sin \alpha x_j\right)}^2\right]\end{array}$$

将 $\alpha =36.795°$ ，$$\beta =78.169°$$ 带入可得：

$$
\begin{gathered}
R^{i}_{j}(\alpha,\frac\pi2-\beta)= \left[\begin{array}{l l l}{cos(36.795 \degree)cos(78.169 \degree)}&{-sin(36.795 \degree)}&{cos(36.795 \degree)sin(78.169 \degree)}\ {sin(36.795 \degree)cos(78.169 \degree)}&{cos(36.795 \degree)}&{sin(36.795 \degree)sin(78.169 \degree)}\ {-sin(78.169 \degree)}&{0}&{c o s(78.169 \degree)}\ \end{array}\right]

\end{gathered}
$$

所以对问题一中理想抛物面公式进行替换，即可得到第二问理想抛物面公式：

$$\begin{array}{r}-\sin 78.169°x_j+\cos 78.169°z_j=F-R-\frac{1}{4b}+b\left[{\left(\cos 36.795°bs78.169°x_j-\sin 36.795°y_j+\cos 36.795°\sin 78.169°x_j\right)}^2+{\left(\sin 36.795°\cos 78.169°x_j+\cos 36.795°y_j+\sin 36.795°\sin 78.169°x_j\right)}^2\right]\end{array}$$

5.3 检验理想抛物面

接下来检验理想抛物面是否会超出标准球面：

图9.垂直角度

图10. $\beta$角度发生偏移

经上图分析可得，理想抛物面是否会超出基准球面仅与$SC$的仰角及角度$\beta$有关，于是本题仅讨论$SC$与$SC$在$XOY$平面上的投影构成的切面上的情况并且角度$\beta$在这个平面上体现$SC$与$X$的夹角。设抛物线上下端点分别为，$F(x_1,y_1)$，$G(x_2,y_2)$，经过坐标转换，可得到图二$F$的坐标为$(x_1^{\prime },y_1^{\prime })$，G点的坐标为$(x_2^{\prime },y_2^{\prime })$。设标准圆半径为$R$，口径为$r$，抛物线口径为$r_2$，那么当$y_1^{\prime }>\sqrt{R^2-r^2}$时，超出标准圆，所以考虑极端情况$y_1^{\prime }=\sqrt{R^2-r^2}$；

此时可以发现，$\angle \beta +\angle \theta +\angle \gamma =\pi$，其中$\angle \theta$是由$R,F$确定的定常数，因此，$\angle \beta$仅与$\angle \gamma$有关，而$\angle \gamma$的变化量正是未旋转时的$CF$与$CI$的夹角；易得向量$\overrightarrow{FC}=\left(-r_a,F-R-\frac{1}{4b}+b{r_a}^2\right)$，向量；$\overrightarrow{IC}=\left(-r,-\sqrt{(R^2-r^2)}\right)$ 于是角度$\mu =arccos\left(\frac{\overrightarrow{FC}\cdot \overrightarrow{IC}}{\lceil \overrightarrow{FC}\rceil \Vert \overrightarrow{LC}\rceil }\right)$ ，即$arccos\left(\frac{r_{3}r-\left(\bar{F}-R-\frac{1}{4b}+b{r}^2\right)\sqrt{\left({\bar{R}}^2-r^2\right)}}{R\sqrt{r_3^2+{\left(\bar{F}-R-\frac{1}{4b}+b{r}^2\right)}^2}}\right)$，所以当且仅当

$$\begin{array}{c}\frac{\pi }{2}-arccos\left(\frac{r_{a}r-\left(F-R-\frac{1}{4D}+b{r}^2\right)\sqrt{(R^2-r^2)}}{R\sqrt{r_a^2+{\left(F-R-\frac{1}{4D}+b{r}^2\right)}^2}}\right)\le \beta \le \frac{\pi }{2}+arcos\left(\frac{r_{3}r-\left(F-R-\frac{1}{4b}+b{r}^2\right)\sqrt{R^2-r^2}}{R\sqrt{r_3^2+{\left(F-R-\frac{1}{4b}+b{r}^2\right)}^2}}\right)\end{array}$$

时，不会超出标准球面。

代入数据得到，当且仅当$1.2568\le \beta \le 1.8848$(弧度）时不会超出标准球面，经计算得出题目二中$\beta =1.364306422991447$，在范围内，符合条件，故不会超出标准球面。

5.4 反射面板调节模型建立

对于反射面板上的主索节点，我们可以通过第一问中的方式，将其沿径向调整至标准曲面上，在通过MATLAB计算后，如图11为模拟工作抛物面，每条边相比于调整之前的边的长度变化比率大小，变化比率越大，边的颜色越蓝，并且发现主索节点间，距离的变化比率最高达到了0.2336695429024%，显然超出了实际中所允许的0.07%，为此，我们需要建立反射面板调节模型。

图11.反射面板调节模型

5.4.1 迭代法

对于这个问题，我们可以使用迭代法，从下至上逐片元调整反射面板，并使其反射的有效信号最大(尽可能多的反射到圆盘内)；对于每一个反射面板的调节，又可以优化为非线性规划问题。

从下至上迭代，当调整点与已经调整好的点不满足约束条件$(>0.07\%)$时，

$$r_{i+1}=r_i+0.01(r_0-r_i)$$

$r_0$表示为调节前的极径，$r_i$表示第$i$次迭代时的极径，其中初值条件$r_0$为理想抛物面下的极径，0.01为迭代步长（因为实际计算中需要尽量满足约束条件$(<0.07\%)$），直至满足条件。

5.4.2 非线性规划

由于面板为三角形，于是光线经面板反射，在$XOY$平面上的投影也是三角形；并且，经计算，投影三角形面积的最小值为$45.473195976021607m^2$，明显大于圆盘面积。因此，使投影三角形和圆尽可能重合，可使反射面板反射到圆盘上的信号量尽量多；

假设原面板坐标为$A_1(x_1,y_{1,}{z}_1)\cdot B(x_2,y_{2,}{z}_2)\cdot C(x_3,y_{3,}{z}_3)$，调整后反射面板坐标为$A_1^{\prime }(x_1^{\prime },y_{1,}^{\prime }{z}_1^{\prime })\cdot B^{\prime }(x_2^{\prime },y_{2,}^{\prime }{z}_2^{\prime })\cdot C^{\prime }(x_3^{\prime },y_{3,}^{\prime }{z}_3^{\prime })$，入射光线，可建立非线性规划如下：

$$\begin{gathered} max\text{反射到圆盘的信号量约束条件} \\ s.t.\frac{||\overrightarrow{A^{\prime }{B}^{\prime }}|-|\overrightarrow{AB}||}{|\overrightarrow{AB}|},\frac{||\overrightarrow{A^{\prime }{C}^{\prime }}|-|\overrightarrow{AC}||}{|\overrightarrow{AC}|},\frac{||\overrightarrow{B^{\prime }{C}^{\prime }}|-|\overrightarrow{BC}||}{|\overrightarrow{BC}|}\le 0.07 \\ |A{A}^{\prime }|,|B{B}^{\prime }|,|C{C}^{\prime }|\le 0.6 \end{gathered}$$

对于一个反射面板，存在以下几种情形：

表 1

情况	待求变量
三个主索节点都没有被固定	三个主索节点都是待求变量
一个主索节点固定，两个主索节点不固定	两个主索节点是待求变量
两个主索节点固定，一个主索节点不固定	一个主索节点是待求变量
三个主索节点都固定	最后检验是否满足约束条件

5.4.3 策略三

是在问题一的基础上增加如下约束条件

$$\frac{\sqrt{(\partial {x_2}^2)+(\partial {y_2}^2)}-\sqrt{(\partial {x_1}^2)+(\partial {y_1}^2)}}{\sqrt{(\partial {x_1}^2)+(\partial {y_1}^2)}}\le 0.07\%$$

即
$$\frac{\sqrt{(\delta x_2{}^2)+(\delta y_2{}^2)}}{\sqrt{(\delta x_1{}^2)+(\delta y_1{}^2)}}\le 1.0007$$

5.5 模型求解

现在我们对这三种方法进行求解，经MATLAB计算，结果如表2，非线性规划比迭代法计算调节后馈源舱接收比更大，从而达到了优化的效果。

表2

方法	反射面积	接收面积	接收比
逐项迭代	66737.9321197487	4040.08026118418	0.060436327885473
非线性规划	66848.5396538344	4040.08049975823	0.060536495085120

最后将方法二非线性规划调节反射面板得到的理想抛物面的顶点坐标，以及调节后反射面300米口径内的主索节点编号、位置坐标、各促动器的伸缩量等结果保存在附录"1.result.xlsx"文件中。

六、问题3模型的分析、建立与求解

6.1 模型分析与建立

模拟平行电磁波射到反射面板发生反射，反射到馈源舱。示意图如图12所示：

图12.反射模型

$\overrightarrow{a}$为平行电磁波，$M_i$为第$i$个主索节点，$\overrightarrow{b}$为反射波，$\overrightarrow{c}=\frac{\overrightarrow{a}}{\overrightarrow{a}}+\frac{\overrightarrow{b}}{\overrightarrow{b}}$为反射面法向量。

以$c$为原点，$\overrightarrow{C{S}_2}$为$y$轴，建立坐标系，得$p(0,-0.534R)$ ，因为$|\overrightarrow{AB}|=1m$，设$\overrightarrow{AB}\perp \overrightarrow{CP}$，得$A(0.5,-0.534R),B(-0.5,-0.534R)$

平行电磁波$\overrightarrow{a}$射到$M_i$形成反射波$\overrightarrow{b_i}$，假设$\overrightarrow{b_i}$经过$A$点，则$\overrightarrow{b_i}=\overrightarrow{M_{i}A}$，取$\overrightarrow{a_0}=\frac{\overrightarrow{a}}{|\overrightarrow{a}|}$，单位向量$\overrightarrow{b_{i_0}}=\frac{\overrightarrow{b}}{|\overrightarrow{b}|}$

可求出$$\overrightarrow{c_i}=\overrightarrow{a_0}+\overrightarrow{b_{i_0}}$$即为第$i$个主索节点的法向量。

可以看出利用$\overrightarrow{a},\overrightarrow{b},\overrightarrow{c}$三向量共面原理与$\overrightarrow{a}$和$\overrightarrow{b}$夹角与$\overrightarrow{b}$和$\overrightarrow{c}$夹角相等，求解出平行电磁波射到反射面板发生反射，反射到馈源舱所在平面的点的坐标。

已知$M_1,M_2,M_3$为相邻的主索节点，因为已求出理想抛物面方程，则已知$M_1,M_2,M_3$调整后的坐标，该问题简化为为平行电磁波被$△M_1{M}_2{M}_3$反射后得$S_{△M_1^{\prime }{M}_2^{\prime }{M}_3^{\prime }}$与圆盘重合部分占反射后形成的$S_{△M_1^{\prime }{M}_2^{\prime }{M}_3^{\prime }}$的比例。示意图如图13所示：

图13

解只需求出$△$三个顶点（主索节点）在馈源舱所在平面的坐标

图14

如图14所示。

设$ M_{1}(x_{1},y_{1},z_{1}),M_{2}(x_{2},y_{2},z_{2}),M_{3}(x_{3},y_{3},z_{3})$ ，$△M_1{M}_2{M}_3$ 的法向量 n ⃗ = M 1 M 2 → × M 1 M 3 → = { a , b , c } \vec{n}=\overrightarrow{M_{1}M_{2}}\times\overrightarrow{M_{1}M_{3}}=\{a,b,c\} n =M1​M2​ ​×M1​M3​ ​={a,b,c}，设$M_1^{\prime }$为$M_1$在馈源舱所在平面的反射点$M_1^{\prime }=(x,y,-0.534R)$，取$\overrightarrow{j}=(0,0,1),\overrightarrow{M_1{M}_1’}(x-x_1,y-y_1,-0.534R-z_1))$

因为三个向量$\overrightarrow{j},\overrightarrow{n},\overrightarrow{M_1{M}_1^{\prime }}$共面，得

$$\left|\begin{array}{ccc}0& 0& 1\\ a& b& c\\ x-x_1& y-y_1& -0.534-z_1\end{array}\right|=0$$

即$a(y-y_1)+b(x-x_1)=0$，由夹角相等，得

$$\begin{gathered} \cos ⟨\vec{j},\vec{n}⟩=\cos ⟨\vec{n},{\overrightarrow{M_1{M}_1}}^{\prime }⟩ \\ \cos \theta =\frac{\vec{j}\cdot \vec{n}}{\left|\vec{j}\right|\vec{n}|}=\frac{\vec{n}\cdot \overrightarrow{M_1{M}_1^{{}^{\prime }}}}{\left|\vec{n}\right|\left|M_1{M}_1^{{}^{\prime }}\right|} \\ c\cdot \sqrt{{\left(x-x_1\right)}^2+{\left(y-y_1\right)}^2+{\left(0.534R+z_i\right)}^2}=a\left(x-x_i\right)+b\left(y-y_i\right)+c\left(-0.534R-z_i\right) \end{gathered}$$
由方程联立求得$x,y$的值，得$M_1^{\prime }$坐标；同理得出$M_2^{\prime }$坐标，$M_3^{\prime }$坐标。

在馈源舱所在平面，求$S_{△M_1^{\prime }{M}_2^{\prime }{M}_3^{\prime }}$与$S_{\text{圆盘}}$重合面积$S_1$

$$\begin{array}{c}{S}_{\Delta M_1^{{}^{\prime }}{M}_2^{{}^{\prime }}{M}_3^{{}^{\prime }}}=\frac{1}{2}\left|\overline{M_1^{{}^{\prime }}{M}_2^{{}^{\prime }}}\times \overline{M_1^{{}^{\prime }}{M}_3^{{}^{\prime }}}\right|\\ S_{\text{圆盘}}=\pi r^2(r=0.5)\\ S_1={\int }_0^{2\pi }d\theta {\int }_0^{1}f(r,\theta )rdr\\ f(r,\theta )=\{\begin{array}{l}1,r\le 1\cap (r,\theta )\in \Delta M_1^{{}^{\prime }}{M}_2^{{}^{\prime }}{M}_3^{{}^{\prime }}\\ 0,\text{其他}\end{array}\end{array}$$

然后推广到每块反射面板用前面建立的单块反射面板调节模型进而求出$k_1=\frac{{\sum }_{i=1}^n{s}_i}{{\sum }_{i=1}^n{\Delta }_i}$，其中$n$为反射面300米口径有效照明区域内反射面板总块数。

$k_1$即为馈源舱的接收比。最后再与$k_2$基准反射球面的接收比相比较$k_z=\frac{{\sum }_{i=1}^N{s}_i}{{\sum }_{i=1}^N{\Delta }_i}$，其中$N$为反射面500米口径区域内反射面板总块数 。

6.2 模型求解

表 3

情况	反射面积	接收面积	接收比
非线性规划	66848.5396538344	4040.08049975823	0.060536495085120
基准球面	193529.689844112	13508.8405906868	0.069802419471494
300米口径有效照明区域	66865.6946103199	4040.08033276622	0.060420823507645

计算得出结果如表3所示，抛物面内的接收比反而比基准球面的接收比小。我们猜想这可能是因为工作抛物面外反射面板的投影三角形面积比工作抛物面内反射面板的投影三角形面积小，因为重叠面积和投影三角形面积不是一个数量级，所以两者近似相等，即两者接收比的分子是相似的，抛物面外的分母却小，所以接收比反而大。

表 4

情况	接收面积
非线性规划	5.880757294300117
基准球面	3.140860402391733
300米口径有效照明区域	5.880757398495225

所以我们又计算了工作抛物面，基准球面以及300米口径有效照明区域的平均每块反射面板的接收面积，如表4所示，发现虽然接收比下降了，但是平均每块反射面板的接收面积增加了。

七、模型评价与推广

7.1 模型评价

本文所建立的"FAST"主动反射面的形状调节模型，是基于观测天体不同方位发射平行电磁波研究，求解理想抛物面进行反射，再利用欧拉法对理想抛物面进行转换，之后先利用迭代法对主索节点进行调整，进而运用非线性规划进行优化，尽量贴近理想抛物面，并求得调节后反射面300米口径内的主索节点位置坐标等结果。最后利用三向量共面原理以及入射角与反射角相等，求解出馈源舱的接收比，并与基准反射球面的接收比作比较。

由于时间精力有限,建立的模型也有一些不足之处，未能实现问题二的方法三的优化，以及在模型建立过程中进行了一些假设，模型考虑的因素有限，求解过程中采取近似计算，对模型的精确性带来一定的影响。

7.2 模型推广

基于观测天体不同方位的"FAST"主动反射面的形状调节策略，求解理想抛物面，利用欧拉法对理想抛物面进行转换，需要注意的一点是如果天体角度过大则会超过基准球面，就需进一步讨论，根据不同的方位对主索节点采用迭代法进行调整，本文研究适用于大部分方位的天体的观测，这就很好的体现了本次建模的意义所在。

第一问中的利用球面与抛物面的各向同性，可以运用于雷达和抗辐射探测器等领域

第二问中的欧拉角法可以应用于经典力学中的刚体研究，与量子力学中的角动量研究。

第三问中的入射角与反射角关系在物理学上，尤其是光学，被广泛应用，本文也是如此。

八、参考文献

1. 钱宏亮.FAST主动反射面支承结构理论与试验研究[D]哈尔滨工业大学2007,4:16-20

2. 李明辉，朱立春.FAST瞬时抛物面变形策略优化分析,贵州大学学报；自然科学版，2012,12:24-25

3. 张帆,曹喜滨,皱经湘.一种新的全角度四元数与欧拉角的转换算法[J],南京理工大学学报；自然科学版,2002,26(4):377-378.

附录

0.支撑材料文件列表

序号	功能	文件名
1	结果计算表格	result.xlsx
2	程序	2.area_获取三角形面积.m
3	程序	3.getarea_获取三角形与圆重叠部分面积.m
4	程序	4.getpoint_获取顶点的投影点.m
5	程序	5.inTra_判断点是否在三角形中.m
6	程序	6.MACM2021_q1_计算问题一中理想抛物线的参数b.m
7	程序	7.MACM2021_q2_计算问题二中所有主索节点位置、促动器伸缩量.m
8	程序	8.MACM2021_q3_计算问题三接收比.m
9	程序	9.mycart2pol_笛卡尔坐标系到球坐标的转换.m
10	程序	10.mypol2cart_球坐标系到笛卡尔坐标系的转换.m
11	程序	11.q1_con_问题一中非线性规划约束条件.m
12	程序	12.q1_fun_问题一中非线性规划目标函数.m
13	程序	13.q1_r1_r2_计算抛物线与圆的极径之差.m
14	程序	14.q2_con_问题二中非线性规划约束条件.m
15	程序	15.q2_fun_问题二中非线性规划目标函数.m
16	程序	16.q2_length_计算单个结点主索节点位置、促动器伸缩量.m
17	程序	17.q2_pro_计算单个反射面的边变化比率.m
18	程序	18.q2_rotate_计算β范围与旋转后的抛物面.m
19	程序	19.q2_shiftPoint1_通过极坐标转换计算根据理想抛物面调整的主索节点坐标.m
20	程序	20.q2_shiftPoint2_通过解方程组计算根据理想抛物面调整的主索节点坐标.m
21	程序	21.q2_shiftPoint3_迭代法调整主索节点.m
22	程序	22.q2_shiftPoint4_非线性规划调整主索节点.m
23	程序	23.q3_fun_问题三中求顶点方程组.m
24	程序	24.q3_receptionRatio_获取问题三的接收比.m
25	参考文献	25.FAST主动反射面支承结构理论与试验研究【1】.pdf
26	参考文献	26.FAST瞬时抛物面变形策略优化分析【2】.caj
27	参考文献	27.一种新的全角度四元数与欧拉角的转换算法【3】.caj

总结

第三问的结果很不理想，后面讨论时发现，我们其实假设了每块反射面板是一个平面，但实际上每块反射面板也是有弧度的，因此导致结果不理想。

相关链接：

• 代码：基于FAST主动反射面的形状调节研究: 2021数模国赛A题 江苏省二等奖 (gitee.com)
  |
  | 24 | 程序 | 24.q3_receptionRatio_获取问题三的接收比.m |
  | 25 | 参考文献 | 25.FAST主动反射面支承结构理论与试验研究【1】.pdf |
  | 26 | 参考文献 | 26.FAST瞬时抛物面变形策略优化分析【2】.caj |
  | 27 | 参考文献 | 27.一种新的全角度四元数与欧拉角的转换算法【3】.caj |

总结

第三问的结果很不理想，后面讨论时发现，我们其实假设了每块反射面板是一个平面，但实际上每块反射面板也是有弧度的，因此导致结果不理想。

相关链接：

• 代码：基于FAST主动反射面的形状调节研究: 2021数模国赛A题 江苏省二等奖 (gitee.com)