How to model (solar) spherical geometry for PV generation?

Modeling
1

Hi there,

I’m using the quap routine from the rethinking package. I’m modeling solar PV generation for a largish region, where there are a few unknown parameters that are easy-ish to constrain.

PV generation can be scaled by cos ( angle of incidence on PV panel ) – the parameter I’m modeling is cos(angle_incidence). This angle is a function of the day of the year, the hour of the day, the solar azimuth (gamma – oriented east or west? zero is due south), the latitude, and the tilt of the PV panel (0 is flat, 90 is upright).

For the raw generation data (actual production), I know the general latitude and the datetime.

My confusion is this: I’m not entirely clear on how to set this up in Stan.

I’ve done this, using quap:

m1 <- quap(
  alist(
    PV ~ dnorm( mu, sigma ),
    lat ~ dnorm(N,39,.5),
    gamma ~ dnorm(N,0,30),
    tilt ~ dunif(N,0,lat),
    sigma ~ dexp(1),
    capacity <- ifelse(solar_altitude(lat,datetime)<0,0,cap),
    cos_incidence <- ifelse(cos_theta_i(solar_altitude(lat,datetime),tilt,gamma,datetime)<0,
                            0,cos_theta_i(solar_altitude(lat,datetime),tilt,gamma,datetime)),
    mu <- capacity*cos_incidence
  ), data = cleandf)

I’ve already used this model to create a prior predictive simulation

Screen Shot 2021-05-09 at 8.17.33 PM
Screen Shot 2021-05-09 at 8.17.33 PM801×637 64.7 KB

Screen Shot 2021-05-09 at 8.19.43 PM
Screen Shot 2021-05-09 at 8.19.43 PM801×637 138 KB

These plots show two things – there are bad data, and the simulated data generally describe the actual data.

I’d like to use Stan to identify the most likely values of the parameters in quap, but when I run it as above, I get this error:
Error in rfelse(1, c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, : could not find function "rfelse"

Am I doing this all wrong? How should models like this be set up?

2

If know neither rethinking, quap or R, but there seems to be a problem with the ternary operator (ifelse). Apart from the data not fitting well anymore, what happens if you remove this thresholding? The error message does sound kind of weird, is that a typo on your part (rfelse) ? Maybe the typo is actually in the R package? The error msg does not sound like something that comes from Stan, so I guess something is going wrong in R?

Besides all this, I think introducing a discontinuity like this is bound to introduce problems.

1 Like
3

Yeah, that’s a good point. I’ll think about how to rewrite it without the ifelse.

4

A quick way forward here might be to use the ulam function from the rethinking package, since ulam uses stan under the hood, but the model is coded in a way similar to quap. You’ll need to make two changes to your workflow to use ulam instead of quap. First, any data transformations (like you’re doing in the ifelse statements), need to be done before passing the data to ulam. Second, the data passed to ulam should include only the columns used in the model.

Another option that’s easier than directly coding the model in stan would be to run your model using the rstanarm package or the brms package. Both of these packages also use stan under the hood, but you code the model using a formula interface, similar to the way glm and other R modeling functions work.

If you want to try brms, Solomon Kurz has translated all of the models in Statistical Rethinking (which uses quap and ulam from the rethinking package) into brms code.

Regarding the error, as Funko_Unko surmised, that looks like an error generated from R due to a typo (rfelse instead of ifelse).

1 Like
5

Actually, I was curious, and without knowing R at all, it looks like the way quap samples from the prior assumes that replacing the first character of whatever comes to the right of the arrow with an “r” lets you sample from it in R.

See rmcelreath/rethinking source: R/map-quap.r

    # function to sample from prior specified with density function
    sample_from_prior <- function( f ) {
        RHS <- f[[3]]
        the_density <- as.character( RHS[[1]] )
        the_rdensity <- the_density
        substr( the_rdensity , 1 , 1 ) <- "r"
        pars <- vector(mode="list",length=length(RHS))
#...
6

I haven’t used quap, but maybe it isn’t set up to handle data transformations within alist and isn’t parsing the ifelse properly.

2 Likes
7

Yes, that’s what I would think as well. Though, data transformations seem to work (search quap function - RDocumentation for ifelse), but probably not arbitrary parameter transformations. The documentation also states (emphasis mine): “Linear models can be specified as formulas …”.

8

In the documentation example, the ifelse is used to define the list of data in the data argument. This is evaluated immediately and doesn’t cause an error. The code in the OP uses ifelse inside alist which delays evaluation and may be causing a problem later when all the expressions in alist are evaluated.

1 Like
9

Ah, I read “inside a list” and didn’t even know that there is such a thing as an alist.

Anyhow, it looks like quap won’t be able to handle this model, and one of the alternatives you provided will be used (y).

1 Like
10

alist doesn’t evaluate its arguments immediately, so they can be processed later, whereas list immediately evaluates its arguments and will throw an error if one of the arguments doesn’t exist or isn’t a valid R language object.

1 Like