Pictured above, you can see that my household has a little solar panel that we use to offset some of the power we need from the utility company. The panel is connected to an inverter which is plugged into a TP-Link smart plug in an outdoor power outlet. The TP-Link smart plug was designed to monitor energy consumption through a proprietary app, but some clever people found out how to reverse engineer the protocol[undefined] and provide a Python client[undefined][undefined] for hackers to use. This allowed us to use the TP-Link device and the Python client to send the incoming power data to an internal database in our house.

To get the biggest bang for our buck, we are interested in orienting the panel for maximal energy production. While we can experiment daily with the different angles and compare different energy production with the database, we may as well start with a simpler model-based method. This post serves to explain just that.

Power

For our initial model, we will assume the radiant flux through the panel only depends on orientation; for now, the different wavelengths of light throughout the day will have no effect on power. To this end, we simply denote power as an inner product $-\langle r, n \rangle$ for coordinate vectors $r, n \in \bbR^3$ associated to the solar radiation and panel's normal vector, respectively. Consider the little diagram below for a silly interactive demonstration.

The radiation vector $r$ should in fact be an expression $r(t, \phi, \theta)$ depending on the time $t$ and location $(\phi, \theta)$ on earth, which we will take to be in the (latitude, longitude) coordinate system. As hinted in a previous post, a nice way to calculate this vector is through the use of frames. Specifically, we introduce a sequence of celestial frame transformations.

With this framework, the center of the solar system is the origin of $F_s$, $o_{F_s} \in \calL$, and a panel on earth at time $t$ and location $(\phi, \theta)$ is the origin of $F_p(t, \phi, \theta)$, $o_{F_p(t, \phi, \theta)} \in \calL$. Thus, the solar radiation is given by the displacement between these two. $$ s(t, \phi, \theta) = o_{F_p(t, \phi, \theta)} - o_{F_s} $$ We calculate the coordinate vector of this in the local earth frame like so. $$ [s(t, \phi, \theta)]_{F_p(t, \phi, \theta)} = -\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} \Big( T_y(t) \cdot T_a \cdot T_d(t) \cdot T_p(\phi, \theta) \Big)^{-1} \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \end{pmatrix} $$

$$ r(t, \phi, \theta) = D^2 \times \frac{[s(t, \phi, \theta)]_{F_p(t, \phi, \theta)}}{\lVert [s(t, \phi, \theta)]_{F_p(t, \phi, \theta)} \rVert^3} $$

Note. The length of $r(t, \phi, \theta)$ is always roughly $1$.

Similarly, the orientation vector $n$ should in fact be an expression of the form $n(\beta, \gamma)$, where $\beta, \gamma \in [-\pi/2, \pi/2]$ are tilting parameters which orient our solar panel. To get these, we introduce a transformation matrix $T_A(\beta, \gamma)$ which rotates $\beta$ along the second coordinate vector and $\gamma$ along the first. $$ T_A(\beta, \gamma) = \begin{pmatrix} \cos\beta & -\sin\beta\sin\gamma & -\sin\beta\cos\gamma & 0 \\ 0 & \cos\gamma & -\sin\gamma & 0 \\ \sin\beta & \cos\beta\sin\gamma & \cos\beta\cos\gamma & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} $$ The frame transformation $F_p(t, \phi, \theta) \rightarrow F_A(t, \phi, \theta, \beta, \gamma)$ associated to this transformation, $$\langle \calL \cup \calD \rangle_{F_p(t, \phi, \theta)} \xleftarrow{T_A(\beta, \gamma)} \langle \calL \cup \calD \rangle_{F_A(t, \phi, \theta, \beta, \gamma)},$$ is such that the third standard basis vector $d_{F_A(t, \phi, \theta, \beta, \gamma), 3}$ would be the panel orientation as in the diagram below. In local earth coordinates, this is calculated as follows. $$\begin{aligned} n(\beta, \gamma) &= [d_{F_A(t, \phi, \theta, \beta, \gamma), 3}]_{F_p(t, \phi, \theta)} \\ &= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} \langle d_{F_A(t, \phi, \theta, \beta, \gamma), 3} \rangle_{F_p(t, \phi, \theta)} \\ &= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} T_A(\beta, \gamma) \langle d_{F_A(t, \phi, \theta, \beta, \gamma), 3} \rangle_{F_A(t, \phi, \theta, \beta, \gamma)} \\ &= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} T_A(\beta, \gamma) \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \end{pmatrix} \end{aligned}$$

The power is now our inner product, up to clamping when the sun is occluded by the earth.

Optimization

We now see that, provided some interval of time $[t_0, t_1]$, the energy produced by a panel at fixed location $(\phi, \theta)$ and orientation $(\beta, \gamma)$ over this time will be given by the following integral. $$ \operatorname{energy}(t_0, t_1, \phi, \theta, \beta, \gamma) = \int_{t_0}^{t_1} \operatorname{power}(t, \phi, \theta, \beta, \gamma) {\rm d}t $$ We would like to find the optimal orientation $(\beta^*, \gamma^*)$. $$ \operatorname{energy}(t_0, t_1, \phi, \theta, \beta^*, \gamma^*) = \max_{-\pi/2\leq\beta, \gamma \leq\pi/2} \operatorname{energy}(t_0, t_1, \phi, \theta, \beta, \gamma) $$

Surely, our model is too complicated to get a closed form expression for the energy integral above. However, we may still implement a quadrature and numerical optimizer to approximate the optimal angle. Below, we implement a Julia[undefined] program which utilizes the Optim.jl[undefined] and QuadGK.jl[undefined] libraries to do just that.

import LinearAlgebra: norm, dot
import QuadGK: quadgk
import Optim: optimize, GradientDescent, Options

# constants
DISTANCE = 1.496e8
RADIUS = 6371
DAYS = 364.24
ALPHA = 23.45 * pi / 180

# timecycle
ω_d(t) = 2 * pi * t
ω_y(t) = ω_d(t) / DAYS

# spherical map coordinate transformation
σ(ϕ) = pi * (ϕ - 90) / 180
η(θ) = pi * (180 - θ) / 180

# frame transformations
T_y(t) = [
  1 0 0 DISTANCE * cos(ω_y(t)) ;
  0 1 0 0;
  0 0 1 -DISTANCE * sin(ω_y(t))
  0 0 0 1
]

T_a = [
  cos(ALPHA) -sin(ALPHA) 0 0;
  sin(ALPHA) cos(ALPHA)  0 0;
  0          0           1 0;
  0          0           0 1
]

T_d(t) = [
  cos(ω_d(t))  0 sin(ω_d(t)) 0;
  0            1 0           0;
  -sin(ω_d(t)) 0 cos(ω_d(t)) 0;
  0            0 0           1
]

T_p(ϕ, θ) = begin
  lat = σ(ϕ)
  long = η(θ)
  [
    -sin(long) cos(lat)*cos(long) sin(lat)*cos(long) RADIUS*sin(lat)*cos(long);
    0          -sin(lat)          cos(lat)           RADIUS*cos(lat)          ;
    cos(long)  cos(lat)*sin(long) sin(lat)*sin(long) RADIUS*sin(lat)*sin(long);
    0          0                  0                  1
  ]
end

T_A(β, γ) = [
  cos(β) -sin(β)*sin(γ) -sin(β)*cos(γ) 0;
  0      cos(γ)         -sin(γ)        0;
  sin(β) cos(β)*sin(γ)  cos(β)*cos(γ)  0;
  0      0              0              1
]

# important vectors in the local earth frame
radiation_m(t, ϕ, θ) = begin
  s = (inv(T_y(t) * T_a * T_d(t) * T_p(ϕ, θ)) * [0 ; 0 ; 0; 1])[1:3]
  -DISTANCE^2 * s / norm(s)^3
end

panel_normal_m(β, γ) = (T_A(β, γ) * [0; 0; 1; 0])[1:3]

power(t, ϕ, θ, β, γ) = begin
  r = radiation_m(t, ϕ, θ)
  r[3] <= 0 ? max(-dot(r, panel_normal_m(β, γ)), 0) : 0
end

# quadrature approximation of energy function
function energy(t0, t1, ϕ, θ, β, γ) 
  quadgk(t -> power(t, ϕ, θ, β, γ), t0, t1)[1]
end

# optimization
function optimize_angle(t0, t1, ϕ, θ; init=[0.0, 0.0], kwargs...)
  objective = action -> -energy(t0, t1, ϕ, θ, action...)
  optimize(objective, init, GradientDescent(), autodiff=:forward, Options(;kwargs...))
end

Note how the program provided an angle that seems to produce optimal energy at some arbitrary location $(\phi, \theta) = (36.1378, -115.1619)$ during a week in winter. Note that the allegedly optimal angle is pointed slightly from vertical toward the south horizon, which makes sense in the northern-hemisphere winter.

References

Related posts