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

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

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?

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

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

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

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.

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

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

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 …”.

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

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

`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