問題的引入

衛星上的電源有太陽能電池電源、化學電源和核電源等多種。太陽能是一種取之不盡、用之不竭的能源。因此，目前絕大多數衛星都採用太陽能電池。這種電池利用半導體材料的光電效應將太陽能直接轉換成電能，能工作幾年甚至幾十年。
早期的太陽能電池片是用半導體矽做成的，現在開始廣泛採用砷化鎵材料，因為它的光電轉換效率高，能達到 20%以上。另外，為了進一步增加太陽能電池片數量，中大型衛星多采用由數塊太陽能電池板連線而成的太陽翼。由於太大，火箭裝不下，所以在發射時太陽翼處於摺疊狀態，待衛星與火箭分離後再展開。為了使衛星上的太陽翼總是朝向太陽，以獲取最大的電能，不少衛星上的太陽翼採用了一些先進技術：一是裝有驅動機構，帶著太陽翼轉動；二是用太陽敏感器來捕獲太陽的方位，從而控制驅動機構轉動太陽翼，使太陽光儘可能垂直於太陽翼，為衛星提供充足的能源。
假設有一顆立方體衛星，我們可以在它的各個面安裝太陽能電池板。衛星環繞地球運動，但軌道任意。怎樣安裝太陽能電池板，才能使衛星吸收最大的太陽輻射功率？本文通過Matlab平臺模擬，對這個問題進行了討論。

分析思路

首先獲取一顆衛星的兩行軌道根數（TLE）資料，提取開普勒軌道根數，以及衛星在近焦點座標系中的座標。在此基礎上，我們將以下列步驟推進這個問題：

1. 利用普朗克黑體輻射公式，計算可見光波段（400nm~800nm）太陽的輻射強度總和M、輻射通量∅、太陽的發光強度I；
2. 利用已知引數，得到太陽在近焦點座標系中的座標；
3. 將太陽在近焦點座標系中的座標轉換為地心慣性座標系中的座標；
4. 再將太陽在地心慣性座標系中的座標轉換到衛星的近焦點座標系中的座標；
5. 利用太陽與衛星的幾何關係和已求出的發光強度I，得到太陽在衛星處的照度；
6. 通過已求出的照度，結合太陽光照方向與衛星各面法向的幾何關係，計算衛星各面太陽能電池板實際最大輸出功率；
7. 分析各方向電池板最大輸出功率資料，提出能夠優化太陽能電池功率 輸出的立方體衛星太陽能電池板安裝或工作方案。

太陽光照強度引數的計算

由普朗克黑體輻射公式，黑體輻射的單色輻射強度為
$$M(\lambda ,T_0)=10^{12}\frac{c_1}{{\lambda }^5(e^{c_2/\lambda T_0}-1)}$$
式中λ為波長(um)，$T_0$為黑體的溫度(K)，通常太陽輻射可以認為是溫度為 5900K 的黑體輻射，$c_1=3.742\times 10^{(}-4)W\bullet u{m}^2$為第一輻射常數，$c_2=14388um\bullet K$為第二輻射常數。單色輻射強度$M(\lambda ,T_0)$單位為$W/(m^2\bullet um)$，在${\lambda }_1-{\lambda }_2$波長範圍內太陽的發光強度總和$M$與輻射通量$\varnothing$分別為
$$M={\int }_{{\lambda }_1}^{{\lambda }_2}M(\lambda ,T_0)\mathrm{d}\lambda ={\int }_{{\lambda }_1}^{{\lambda }_2}10^{12}\frac{c_1}{{\lambda }^5(e^{c_2/\lambda T_0}-1)}\mathrm{d}\lambda$$
$$\varnothing =4\pi {R_S}^{2}M$$
式中$R_s=6.9599\times 10^{8}m$為太陽半徑，$M$的單位為$W/(m^2)$。假設太陽所發出的總輻射通量在空間方向上的分佈是均勻的，則在${\lambda }_1-{\lambda }_2$波長範圍內太陽的發光強度為
$$I=\varnothing /4\pi =R_s^{2}M$$
根據距離平方反比定律，在${\lambda }_1-{\lambda }_2$波長範圍內太陽在衛星處的照度(即大氣層外的照度)$E_s$為
$$E_s=I/D^2=R_s^{2}M/D^2$$
根據上式，若知道太陽與衛星的距離，便可得到太陽在衛星處的照度。但根據實際情況，衛星可能執行到地球陰影面，導致光照消失。在後續計算照度時應注意這一點。

太陽在各座標系中座標的轉換

考慮到計算照度時需要得到太陽到衛星的距離，將太陽座標轉換到衛星的近焦點座標系將會便於求解。
太陽軌道引數：

已知以上引數後，便能夠求解太陽在其自身的近焦點座標系中的座標。太陽在近焦點座標系中的座標為：
$$x_{sun,o}=a(cosE-e)$$
$$y_{sun,o}=a\surd (1-e^2)sinE$$
$$z_{sun,o}=0$$
近焦點座標系到地心赤道慣性座標系的轉移矩陣可表示為：

因此太陽在地心赤道慣性座標系中的座標為：

用同樣的方法，我們可以得到衛星的近焦點座標到地心赤道慣性座標系的轉移矩陣${T_{ip}}^{\ast }$。設$Y_{isun}$為太陽在衛星近焦點座標系中的座標，$Y_{p}sun$為太陽在地心赤道慣性座標系中的座標，則
$$Y_{p}sun={T_{ip}}^{\ast }\ast Y_{i}sun$$
利用矩陣左除的方法，可由$Y_{p}sun$和${T_{ip}}^{\ast }$得到$Y_{i}sun$，從而達成太陽座標轉換的目標。

太陽在衛星處照度的計算

在照度的計算中，最關鍵的是確定當前時刻衛星所在位置是否能夠收到光線，即衛星是否處於地球陰影面。
任意時刻，太陽、地球、衛星都處於同一個平面，它們的相對位置關係如下圖：

圖中，藍色的圓表示地球，紅色的圓表示太陽。衛星在圍繞地球的軌道上運動，圖中黑色的點A標出的是衛星有照度-無照度的臨界位置：若衛星到太陽的距離大於該位置到太陽的距離，則衛星進入地球陰影區，無光照。
由於太陽到地球的距離遠遠大於地球半徑，可近似認為太陽光是平行光，這樣，點A的位置可由幾何關係確定。設D是過A向地球圓所作切線與地球圓的焦點，則$AD\Vert OC$。
延長OC，過A向OC作垂線，垂足為B，則AB=OD=地球半徑；設衛星在近焦點座標系中的位置座標為$(x,y,z)$，則$AO=\surd (x^2+y^2+z^2)$；由於ABO是直角三角形，可以根據勾股定理計算出BO的長度；OC表示太陽到地心的距離，可由太陽的座標直接求出，得到BC=BO+OC；同時ABC也是直角三角形，AB、BC長度已知，根據勾股定理便能夠計算出AC的長度，而AC就表示了臨界點衛星到太陽的距離$D_0$。
在衛星執行的每一時刻，都對衛星到太陽的距離D進行計算，並將計算結果與臨界點的距離值進行比較。若$D>D_0$，說明衛星處在陰影區，照度為0；若$D<D_0$，說明衛星處在光照區，依照公式$E_s=I/D^2=R_s^{2}M/D^2$計算照度即可。

衛星各面太陽能電池板實際最大輸出功率的計算

立方體衛星的各面示意如下：

由於衛星在軌道座標系統中俯仰、偏 航、翻滾角均為0度，故衛星六個面中，必有一面法線平行於軌道切向（A），一面法線指向地心（B）。與這兩面方向不同的另一面（C），其法線可由這兩面法線的叉乘得到。剩下的三個面與這三個面相對，法向剛好相反。也就是說，有了A面和B面的法線方向，便能求出全部6個面的方向。
A面的法向，即軌道切向，可由兩相鄰時刻衛星的運動距離向量得到。當時間差選取地足夠小，衛星在兩相鄰時刻運動的距離向量，近似平行於衛星在該位置處的切線；由於計算在衛星的近焦點座標系中進行，B面的法向就是衛星的位矢，由衛星的座標直接確定。
太陽光線方向可由太陽與衛星的座標差確定。設太陽光線向量a,衛星某一面法向為b，a與b的夾角為θ,則

如果一塊太陽能電池板面積為S，電能轉化效率為η，最大輸出功率為$W_{max}$，與太陽的距離為D，那麼太陽能電池板實際最大輸出功率為
$$W=max(0,min(\eta E_{s}cos\theta ,W_{max}))$$

Matlab程式碼講解

1）計算太陽光照強度與位置引數
依照光照強度引數計算及座標轉換中的資料及公式，代入計算即可。
太陽位置引數計算：

asun=149597870;                                 %半長軸
esun=0.016709-0.4204*0.1^3;                     %偏心率
isun=0;                                         %軌道面傾角
Msun=pi*(357.5277233+35999.050)/180;            %真近地點角
omegasun=pi*(282.9404+1.7201)/180;              %近地點幅角
Esun=Msun+esun*(1+esun*cos(Msun))*sin(Msun);    %偏近點角
ibusasun=pi*(23.439291-0.0130042)/180;          %黃赤夾角

400-800nm太陽發光強度計算：

syms lamda
T0=5900;c1=3.742*10^(-4);c2=14388;
f=10^12*c1/(lamda^5*exp(c2/(lamda*T0))-1);
M=double(int(f,lamda,0.4,0.8))*10^(-6);
Rs=6.9599*10^8;
I=Rs^2*M;

需注意對符號函數進行定積分運算時，最好對計算結果進行double強制型別轉換，否則可能因型別不匹配而報錯。

對太陽進行座標轉換
首先根據公式計算太陽在它自己的近焦點座標系中的座標：

xsunp=asun*(cos(Esun)-esun);
ysunp=asun*sqrt(1-esun^2)*sin(Esun);
zsunp=0;

然後利用轉移矩陣，轉換到赤道慣性座標系：

T1=[1 0 0;0 cos(ibusasun) -sin(ibusasun);0 sin(ibusasun) cos(ibusasun)];
T2=[cos(omegasun) -sin(omegasun) 0;sin(omegasun) cos(omegasun) 0;0 0 1];
Tsun=T1*T2;
[Ysuni]=Tsun*[xsunp;ysunp;zsunp];

自行定義了函數ECI2PF，用以完成從地心慣性座標系到衛星的近焦點座標系的座標轉換：

function [Y1] = ECI2PF (Y0, Kep)

% 開普勒引數
a = Kep(1);  Omega = Kep(4);
e = Kep(2);  omega = Kep(5);
i = Kep(3);  M0    = Kep(6);

%轉換矩陣計算
Tip1=[cos(Omega),-sin(Omega),0;sin(Omega),cos(Omega),0;0,0,1];
Tip2=[1,0,0;0,cos(i),-sin(i);0,sin(i),cos(i)];
Tip3=[cos(omega),-sin(omega),0;sin(omega),cos(omega),0;0,0,1];
Tip=Tip1*Tip2*Tip3;

% 座標系轉換
Y1=Tip\Y0;
end

函數的輸入為赤道慣性座標系中的座標Y0和衛星的開普勒引數Kep，輸出所求點在衛星的近焦點座標系中的座標Y1。Tip原本是近焦點座標系-赤道慣性座標系的轉移矩陣，這裡利用MATLAB帶有的「左除」這一特殊運算，完成了反向轉換的功能。
通過該函數，就把太陽的座標從赤道慣性座標系轉到了衛星的近焦點座標系：

[Ysun]=ECI2PF(Ysuni,Kep);

照度的計算
太陽座標已有，現在確定衛星的座標：

[Y0] = State (gm, Kep, dt);
Ys=0.1^3*Y0(:,1);

State函數在上次實驗中已有使用，其功能是確定衛星在它自己的近焦點座標系中的位置和速度座標，在這裡不再詳細闡述。第二行命令的作用是提取衛星的位置座標，並將其存放在向量Ys中。
根據日-地-衛星圖中說明的原理，衛星是否處於地球陰影區取決於它與太陽的距離：

ds_sun0=6371^2+(sqrt(sum(Ysun.^2))+sqrt(sum(Ys.^2)-6371^2))^2;
ds_sun=sum((Ys-Ysun).^2);

ds_sun0表述臨界距離，ds_sun表示衛星與太陽的實際距離。當實際距離大於臨界距離時，照度為0.否則按照照度計算公式計算衛星處的照度：

if ds_sun>ds_sun0
    Es(dt/60)=0;
else
    Es(dt/60)=I/ds_sun;
end

對上述步驟按一定時間步長進行迭代（控制總時長12000s），即可得到從規定時刻起200分鐘內太陽在衛星處的照度曲線：

for dt=60:60:12000
…
end

計算衛星各面段太陽能板實際最大輸出功率
這一步關鍵在於找到3個互相垂直的法線方向與太陽光照方向（太陽與衛星連線方向）的夾角餘弦值。
首先是衛星前進方向，法線方向與軌道切線平行，利用很短時間內衛星的位移方向近似替代。
在前一時刻，把衛星座標存放在臨時變數temp中：

temp=Ys;

到下一時刻，用現在的衛星座標與temp相減，就是短時間內衛星的位移：

Yds=Ys-temp;

利用向量夾角公式，得到該方向的夾角餘弦值：

costheta1=sum(Ysun.*Yds)/sqrt(sum(Ysun.^2)*sum(Yds.^2));

其次是法線指向地心的方向。由於各座標均屬於衛星的近焦點座標系，衛星的座標就表示了該法線方向，運算相對方便：

costheta2=sum(Ysun.*Ys)/sqrt(sum(Ysun.^2)*sum(Ys.^2));

最後一個法向用前兩個法向的叉乘表示：

Yon=cross(Ys,Yds);
costheta3=sum(Ysun.*Yon)/sqrt(sum(Ysun.^2)*sum(Yon.^2));

利用三個餘弦值，結合太陽能電池板實際最大輸出功率的計算公式，便能得到各方向太陽能板實際最大輸出功率（有效面積60平方釐米，效率30%，最大輸出功率2.4W）：

W1(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*costheta1,2.4));
W2(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*(-costheta1),2.4));
W3(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*costheta2,2.4));
W4(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*(-costheta2),2.4));
W5(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*costheta3,2.4));
W6(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*(-costheta3),2.4));

模擬結果分析

1）確定時刻起至之後 200 分鐘內太陽在衛星處的照度曲線
衛星：amateur；起始時刻：2020.10.4.22:10

曲線形如矩形方波。由於衛星軌道的大小範圍遠小於日地距離，故衛星在非陰影區運動時，其到太陽的距離變化對照度的影響可以忽略不計，近似為恆定值；而一旦運動到陰影區，照度立馬減小到0。這說明模擬結果是合理的，與實際情況能夠較好契合。

2）立方體衛星各面太陽能板實際最大輸出功率曲線
衛星各面如下圖所示：

以下為各面最大輸出功率在2020.10.4.22:10到此後200分鐘的曲線圖及各圖的分析：
正面與背面

上圖中有兩個完全相同的脈衝波形，兩波形之間輸出功率恆為零。由於衛星的繞轉具有周期性，可估算amateur衛星的繞轉週期約7000s。
波形在靠近2000s、超出8000s處有兩個陡增點，且這兩個陡增點與照度曲線的陡增點完全匹配，說明衛星在這兩個時刻從陰影區進入了光照區，並受到一定功率的太陽光照射；而後曲線緩慢上升，又逐漸減小到0，這也是與實際情況相符合的：

位置1：衛星進入光照區，照度突變為一定值；
位置2：θ=0，功率達到最大；位置1-位置2 θ減小，功率增大
位置3：θ=90，功率減小至0；位置2-位置3 θ增大，功率減小；位置3-位置1過程由於正面始終背向太陽光，功率保持0不變。

背面與正面的輸出功率曲線具有對稱性：位置3（圖7中大於4000s處的起始點）是正面輸出功率減小到0的點，卻是背面輸出功率從0開始逐漸增加的點；從位置3開始，背面輸出功率逐漸增加到最高點（位置4），又逐漸下降，指導進入地球陰影區（位置5），功率突變減小至0。
右面與左面

背向地球的右面，最大輸出功率曲線始終是平滑的，說明這一面未受到陰影區的影響：從位置2-位置3，右面法線與陽光方向夾角減小，功率增大；位置3-位置4，右面法線與陽光方向夾角增大，功率減小；而在受到陰影區影響的位置2-位置4過程，由於幾何遮擋，功率始終為0。

假設沒有地球的存在，左面的功率曲線應與右面在一個週期內實現對稱。但因為地球的遮擋，左面本應到達最大功率的位置5-位置1區間，完全沒有照度，在影象上體現為在原本的脈衝峰值處呈現矩形凹陷，且凹陷點的時刻與照度曲線相匹配。
上面與下面


對amateur衛星來說，上面最大輸出功率始終為零，下面最大輸出功率則同照度曲線完全一致。若更換不同的衛星，則上面、下面的功率曲線可互換，但呈現的規律保持不變。也就是說，對於這兩個面，在非陰影區，必有一個面始終受到光照且其法線方向同陽光方向的夾角保持不變，而另一個面始終處於被幾何遮擋的狀態。這與實際情況也是相符合的。

總結–立方體衛星太陽能電池板安裝方案的設計

對於立方體衛星太陽能電池板的安裝方案，有以下假設：
① 為每個面都安裝電池板是不現實的；
②電池板的安裝應滿足對稱性，即某一面若安裝電池板，其對應面也應安裝電池板；
② 最優安裝方案應使得所有電池板獲得的總功率在可選方案中最大。
為選定最優安裝方案，我們需要求出各面的最大輸出功率，或者是最大輸出功率的相對大小，並對結果進行排序。這一步的程式實現為：

S1=sum(W1);
S2=sum(W2);
S3=sum(W3);
S4=sum(W4);
S5=sum(W5);
S6=sum(W6);
[S,I]=sort([S1+S2 S3+S4 S5+S6],'descend')

W1-W6表示分別正、背、右、左、上、下面在各時刻的最大輸出功率，則S1-S6就能夠代表各面平均最大輸出功率的相對。由於假設②，進行排序的是正+背、右+左、上+下面平均功率之和。一組面的平均功率之和越大，在其上安裝太陽能電池板的效率越高。
使用不同衛星得到的結果如下：

衛星名稱	正+背面平均功率相對大小	右+左面平均功率相對大小	上+下面平均功率相對大小	排序（降序）
amateur	119.2251	103.6662	190.2461	上+下＞正+背＞右+左
nnss	154.4545	115.8150	103.6619	正+背＞右+左＞上+下
noaa	132.4896	112.3542	143.6318	上+下＞正+背＞右+左
radar	91.2918	85.8307	290.1408	上+下＞正+背＞右+左
geo	68.3337	307.3655	49.9478	右+左＞正+背＞上+下
gps-ops	131.7696	270.6468	72.5789	右+左＞正+背＞上+下

選取多組衛星並進行排序比較後，可以得到這樣的結果：

① 正+背面，即法線平行於軌道切線的兩個面，其平均功率總不會最小；
②右+左與上+下面，其平均功率與具體的衛星有很大的關係：一般來說，高緯度執行的衛星，其上+下面獲得的功率最大；赤道平面衛星或低緯度執行的衛星，其右+左面獲得的功率最大。
據此，一個立方體衛星若有兩組方向可以安裝太陽能電池板，其最佳安裝方案為：
①正+背面必安裝電池板；
②若衛星執行軌道緯度較高，在上+下面安裝電池板；若為赤道平面衛星或低緯度執行的衛星，在右+左面安裝電池板。