先に買ったAliの150Wソーラーパネル。
- 表面がPET樹脂で燃えるなんてひどいぜ
というのはとりあえず措くとしても、
- マジで150Wあるの？
という疑惑と、
- 正面から外れたときの発電量の低下が妙に急では？
という疑問が頭の中を渦巻いている。
こうした状況では、パネルと太陽が正対した状態を100%とおき、そこからの偏角に応じて変換効率がどのように変化するかプロットしておけば役に立ちそう。
もちろん、インバータの不安定さを解決しないといろいろ不正確になるではあるんだけど、*1 大量のデータで統計的に推測すれば、わりによいguessができる気がする。
パネルは固定、太陽はPyEphemが方角+仰角で時間に応じた角度を返してくれるので、これは楽勝。
と思ったんですけどね。甘かったw というか、ナメてた。線形代数を30年以上使ってなかったら完全に忘れてる。
それでもまあ、なんとなくの勘はある。まず、2つの空間的な方角のなす角を求めるには、ベクトルのなす角を使うのが楽そう。
これは2つのベクトルをa、bと置くと、内積a・bの定義:
a・b = |a||b|cosθ
より
cosθ = a・b / |a||b|
となるから。cosθから角度θ（シータ）はアークコサインで出る。
θ = arccos(cosθ)
これで求まるはず。Python的には:
import numpy as np # numpyでベクトルを扱いやすく
def vecs2deg(vec1, vec2):
cos_theta = np.dot(vec1, vec2) / (np.linalg.norm(vec1) * np.linalg.norm(vec2))
return rad2deg(np.arccos(cos_theta))
とかやると、vecs2deg(ベクトル1, ベクトル2) により、ベクトル1とベクトル2がなす角が返ってくる。
と、ここまでは楽勝だったんだけど、角度をベクトルに直すのが実にめんどかった。
太陽の位置とかパネルの方向とかは方角と仰角だけで示されている。原点からの角度と距離で座標を示す極座標系ってのがあるんだけど、方角と仰角だけで示されてる物体というのは、(x, y, z)の直交座標系ではなく、(距離, 角度, 角度)で示される極座標系（球座標）に乗ってると考えた方が抽象化がうまくいく。
球座標系と直交座標系は相互に変換できる。その一般的な変換は、x軸となす角をφ（ファイ）、z軸となす角をθ、原点からの距離をrとすると
x = r sinθ cosφ y = r sinθ sinφ z = r cosθ
である。これをPythonで素直に実装すると:
import numpy as np
from numpy import sin, cos # 三角関数をそのまま使いたい
def polar2vec(phi, theta, r) -> np.array:
x = r * sin( theta ) * cos( phi )
y = r * sin( theta ) * sin( phi )
z = r * cos( theta )
return array([x, y, z])
である。
では南中してる太陽の方角を長さ1のベクトルであらわしてみよう。南中だから方位角は180度。ここ数日の太陽の南中高度（仰角）は58度なので
>>> polar2vec(180, 58, 1) array([-0.59419463, -0.79544254, 0.11918014]) >>>
さっそくヘンな値が出た！ なんでこんなに微小な数字に？？
もちろん、ちょっと調べたらわかる。numpyでは角度にはラジアンを使うのが普通なのだ。
では角度をラジアンに変換して計算するようにしてみよう。
import numpy as np
from numpy import sin, cos
def polar2vec(phi, theta, r) -> np.array:
ph = np.deg2rad(phi)
thet = np.deg2rad(theta)
x = r * sin( thet ) * cos( ph )
y = r * sin( thet ) * sin( ph )
z = r * cos( thet )
return array([x, y, z])
こいつで改めて南中太陽の方角ベクトルを作ってみる。
>>> polar2vec(180, 58, 1) array([-8.48048096e-01, 1.03855939e-16, 5.29919264e-01]) >>>
e-01とかe-16とかは浮動小数点数で、e-01は10-1で1/10、e-16は10-16で1/10000000000000000だ。要するに「ほぼ0」である。yがほぼ0。
φはx軸とのなす角なので、南を180度とするのであれば、
- x軸: 南北方向の軸。数字が増えると北に行く。
- y軸: 東西方向の軸。数字が増えると西に行く。
- z軸: 上下方向の軸。数字が増えると上に行く。
ということになる。つまりこのベクトルは、正面に0.848進み、左右には移動せず、上に0.53進んだ方角に太陽がある、という意味になる。
58度ってほとんど60度だから、正三角形の半分の直角三角形の斜辺くらいになる。ということは正面に1進んで上に√3くらい進んだくらいになるはずだ。√3は1.732くらいだから、前進する成分より上に上がる成分の方が1.7倍くらいになってるはず。数字はそうなっていない。
と考えていくと、太陽高度は仰角、極座標系はz軸となす角である。ここが間違ってる。
一般形を崩さず、可読性を犠牲にせずに仰角をz軸とのなす角に直したいから、ちょう簡単な関数を書こう。仰角elevetion angleを極座標のθに直す関数なので、名前はea2theta。単に(90-仰角)を返せばよい。
def ea2theta(elev_angle):
return 90 - elev_angle
やってみよう。
>>> polar2vec(180, ea2theta(58), 1) array([-5.29919264e-01, 6.48963931e-17, 8.48048096e-01]) >>>
xとzが入れ替わり、0.53進んで0.85上がるようになった。0.848/0.53は1.6。1.732にちょっと遠いけど、だいたいこんなもん？？
気持ちわるいので、もうちょっと検算しよう。
南に向かって1進むだけのベクトルは真南の地平線の角度を示すベクトルだ。これと太陽方向のベクトルを、最初に書いた関数に入れて角度を出してみる。太陽と地平線の角度はすなわち太陽高度だ。
>>> vecs2deg(np.array([-1, 0, 0]), polar2vec(180, ea2theta(58), 1)) 58.00000000000001 >>>
OKだった！！
それでは、この太陽と、庭のパネルなす角を比べよう。庭のパネルはブロック塀に立てかけてあり、40度の傾きにしてある。方角は195度（南南微西）である。
設置角は地面（xy平面）に対する立ち上がりなので、パネル正面方向の角度は設置角0度のとき真上を向き、つまりz軸に等しい。角度が付けばそれはそのままz軸とのなす角となるので、つまりそのまま極座標系のθとして使えるわけだ。
また方位角については、は北から東〜南〜西と回っていくので右回り。極座標系の方位はx軸とのなす角なので左回り。0度方向が共通で絶対値も共通。符号だけが逆になる。
これも簡単な関数にしよう。方位角アジマスを極座標のφに直すので名前はaz2phi。
def az2phi(azimuth):
return -azimuth
確認しよう。
>>> polar2vec(az2phi(191), 40, 1) array([-0.63097779, 0.12264966, 0.76604444]) >>>
xが正の数（南）、yが小さな正の数、zが大きな正の数。よし！
それでは両者のなす角を求めよう。
>>> v_suntran = polar2vec(180, 58, 1) >>> v_panel = polar2vec(az2phi(191), 40, 1) >>> vecs2deg(v_suntran, v_panel) 19.77284232715524 >>>
やっとできた〜。
しかしまだやることがある。PyEphemの返す角度がラジアンだったり度数だったりするんだけど、まだイマイチ慣れてなくて、どっちが返ってくるかわかんないんだよね。
>>> obs.date = datetime.datetime.utcnow() ; sun.compute(obs) # 時間更新 >>> sun.alt # 太陽高度に直接アクセスすると 0.9489879608154297 # ラジアン >>> sun.az # 太陽方位角に直接アクセスすると 2.5738658905029297 # やっぱりラジアン >>> print(sun.alt) # ところがprint()に食わすと 54:22:22.8 # 度数…だけど角度分や秒がめんどい
もうnumpyに入れよう。
>>> np.rad2deg(sun.alt) 54.37300496345047
角度の浮動小数点数。これが欲しかった。
続いて「現在の太陽方向」をベクトルにしてみる。
>>> polar2vec(az2phi(np.rad2deg(sun.az)), ea2theta(np.rad2deg(sun.alt)), 1) array([-0.49112567, -0.31322327, 0.8128264 ])
ちょっと複雑になったけど、まあ可読性は保たれてるからよし。
時間を更新しながらパネルとの角度を見ていこう。
>>> obs.date = datetime.datetime.utcnow() ; sun.compute(obs) >>> v_sun = polar2vec(az2phi(np.rad2deg(sun.az)), ea2theta(np.rad2deg(sun.alt)), 1) >>> vecs2deg(v_panel, v_sun) 25.084119443728053 >>> obs.date = datetime.datetime.utcnow() ; sun.compute(obs) >>> v_sun = polar2vec(az2phi(np.rad2deg(sun.az)), ea2theta(np.rad2deg(sun.alt)), 1) >>> vecs2deg(v_panel, v_sun) 25.01952869385647 >>> obs.date = datetime.datetime.utcnow() ; sun.compute(obs) >>> v_sun = polar2vec(az2phi(np.rad2deg(sun.az)), ea2theta(np.rad2deg(sun.alt)), 1) >>> vecs2deg(v_panel, v_sun) 24.97078892865783 >>>
時々刻々と変わるようになった。
では3秒ごとに更新するようにしよう。
def main():
while True:
obs.date = datetime.datetime.utcnow()
timestamp = ephem.localtime(obs.date).strftime('%Y-%m-%d %H:%M:%S')
sun.compute(obs)
v_sun = polar2vec(az2phi(np.rad2deg(sun.az)), ea2theta(np.rad2deg(sun.alt)), 1)
print('%s, %0.2f' % (timestamp, vecs2deg(v_panel, v_sun)))
time.sleep(3)
>>> main()
2023-03-09 11:54:22, 20.18
2023-03-09 11:54:25, 20.17
2023-03-09 11:54:28, 20.16
2023-03-09 11:54:31, 20.15
2023-03-09 11:54:34, 20.13
2023-03-09 11:54:37, 20.12
2023-03-09 11:54:40, 20.11
2023-03-09 11:54:43, 20.10
2023-03-09 11:54:46, 20.09
2023-03-09 11:54:49, 20.08
2023-03-09 11:54:52, 20.07
ようやくできたっぽい。工学系のマトモな学生なら一瞬でできそうなのに、何時間かかってるやら。
とはいえ出来たものは便利だ！ いまはdatetime.datetime.utcnow()で取った現在時刻を使ってるから時々刻々のズレを記録するようになってるけど、太陽の位置は日付時刻がわかれば特定できるので、これまで取ったデータにも角度情報を付加することができるのである。
とか自画自賛してたら、「オイラー角」という声が聞こえてきた。
あっこれだっ！（泣）
*1:力率すら狂っていくことが昨日のデータでわかっちゃって脱力してる。庭に下ろしてるインバータで、今朝起動せず。パネルの抜き差しで起動…こんなのデプロイしたくないよ。
