Skip to content

Differentiable Physics

Status: durable method explainer. 🚧 Planned / 🔬 Research — none of this is implemented yet. This page explains the method: what differentiable physics (DP) is, the inversion idea at its core, and the two reusable building blocks (a single-site solar model and an aggregate-fleet model). The plans for applying it live in the roadmap: DP is one of the candidate estimators for metered-generator effective capacity (v0.7), and the engine behind the graph-structured disaggregation of net substation power (post-v2 research). The tooling counterpart to this page — which estimation problems belong in convex optimisation (CVXPY) rather than PyTorch, and the cvxpylayers bridge between the two — is Convex Optimisation. The Python in this document is illustrative sketch code, not the implementation. See the roadmap index for status conventions.

Why differentiable physics?

Pure black-box models (like standard neural networks or XGBoost) require vast amounts of data to learn basic physical invariants and are prone to overfitting or hallucinating non-physical behaviour. Differentiable physics (DP) instead injects a mechanistic physical model directly into the computational graph, for three core reasons:

  • Data and sample efficiency: The model does not have to spend capacity re-learning solar geometry (where the sun is) or the shape of a turbine power curve (roughly cubic in wind speed between cut-in and rated, then flat to rated power, then zero beyond cut-out). The physics engine supplies these constraints for free, letting the ML components focus entirely on atmospheric and behavioural nuances.
  • True invertibility (inputs as parameters): Because every operation in the physics engine preserves gradients, we can backpropagate errors all the way to the inputs. This lets us treat unobserved variables (like local irradiance) or system configurations as parameters that can be solved via gradient descent.
  • Interpretability: Instead of inspecting uninterpretable latent layers, the system updates explicit physical parameters like tilt, azimuth, or capacity. This lets engineers immediately audit the model's assumptions.

It is established practice to add a learned residual on top of a physical generator model, and adding the physics makes the forecast interpretable without making it less accurate — the energy-forecasting review reports Gijón et al. (2025)'s turbine result and its limits: the gain is measured from wind at the moment predicted, not from a multi-day forecast, and nobody has put a differentiable generator model inside a distribution network's probabilistic net-demand forecast.

Graceful degradation when an input is missing

A fourth reason, less often stated: a physical forward model has a defined output for any input state. Where an input is absent — a missed NWP (Numerical Weather Prediction) run, an unavailable variable, stalled telemetry — you substitute a climatological prior or a physical bound and let the physics propagate it. There is no branch, no fallback path, and no if data_is_missing:. The same code does the right thing because of how it is arranged.

The behaviour that falls out is exactly the degradation behaviour this project wants from its forecasts. Physics supplies the envelope, the learned residual sharpens it, and as data degrades the residual head has less to work with, so the answer relaxes toward the prior: wider, still bounded, and still true. Contrast a purely learned model, which under an unfamiliar missingness pattern rides default-routing decisions that were never evaluated against any data.

This graceful-degradation behaviour makes missingness robustness a design requirement of any DP estimator we build, rather than a nice-to-have: the priors and bounds substituted for each absent input have to be chosen deliberately, and the estimator's uncertainty must widen honestly when it leans on them. It is scored that way in the capacity-estimation head-to-head. The wider principle is Inherent Stability.

The core idea: inversion through a differentiable forward model

The fundamental insight is to treat a meter reading as the output of a forward model that takes the unobserved quantities as inputs. For a substation:

\[P_{\text{obs}}(t) = P_{\text{demand}}(t) - P_{\text{pv}}(t) - P_{\text{wind}}(t) + P_{\text{loss}}(t)\]

where \(P_{\text{obs}}\) is the observed meter reading, \(P_{\text{demand}}\) the latent (unobserved) demand, and \(P_{\text{pv}}\), \(P_{\text{wind}}\), \(P_{\text{loss}}\) are the photovoltaic (PV) generation, wind generation, and network losses respectively.

Each right-hand-side term is modelled explicitly as a physical function of weather and time, with its own latent parameters (capacities, orientations, power curves). The forward model is differentiable end-to-end: every component is implemented in a differentiable framework (JAX or PyTorch), so that gradients of the reconstruction error with respect to all latent parameters and states can be computed and used for learning. Training minimises the discrepancy between the forward model's output and the observed meter readings, jointly optimising every latent term simultaneously. This is the "inversion" step — running the forward model backwards, constrained by physics.

The pay-off is not just the parameter estimates but the latent terms themselves: subtract the reconstructed generation and what remains is the latent demand, free of the confounding effect of distributed energy resources (DERs). How this project deploys the idea — which terms are metered vs unmetered, which DER types are tractable, and the full graph-structured engine over the substation network — is planned in Net-demand disaggregation.

A note on the name: "differentiable physics" is shorthand for a physics simulator written in an autodiff framework so that it composes with gradient-based learning — the differentiability is the enabling property, not the point. The underlying activity is inverse modelling, and it is shared with the convex route: both toolchains write a forward model and invert it against observations, and both host physics. What actually separates them is the expressiveness of the forward model versus the guarantees on the inversion — unpacked in Two routes to the same inverse problem.

The core building block: DifferentiableSolarPlant

We model each physical parameter as a learnable Normal distribution \(\mathcal{N}(\mu, \sigma^2)\) — a mean-field variational posterior — and train with the reparameterisation trick (rsample()) so gradients flow through the sampling step.

Crucially, the training objective is an ELBO, not a bare reconstruction loss: a power-reconstruction term plus a KL term that pulls each posterior toward a fixed physical prior. ELBO stands for evidence lower bound, and KL for Kullback-Leibler. The KL term is not optional. Minimising power error alone always rewards shrinking \(\sigma \to 0\), so the parameter "uncertainty" we are trying to capture would simply collapse. The prior does double duty: it keeps the posterior spreads honest, and it injects weak domain knowledge (e.g. "panels point roughly south at a typical UK roof pitch") that regularises sites with little data.

Two practical details the ELBO must get right — both classic failure modes of hand-rolled variational inference:

  • The reconstruction term must be a proper likelihood, not a bare MSE. Bare MSE silently assumes an observation-noise scale of 1 in whatever units the power happens to be in. This assumption makes the balance between reconstruction and KL arbitrary — the posterior then either collapses onto the prior or ignores it, depending on nothing but the units. Use a Gaussian likelihood with a learnable observation-noise scale \(\sigma_{\text{obs}}\) (see negative_elbo in the sketch).
  • Scale the KL for minibatches. The KL term regularises the whole-dataset objective once, so on a minibatch it must be weighted by batch_size / dataset_size — otherwise the effective prior strength depends on the batching.

Here's a quick Python code sketch of the rough idea for applying differentiable physics to solar:

import torch
import torch.nn as nn
import torch.distributions as dist

STC_IRRADIANCE = 1000.0  # reference plane-of-array irradiance at standard test conditions (W/m^2)
GROUND_ALBEDO = 0.2  # typical broadband ground reflectance
STC_CELL_TEMP = 25.0  # cell temperature at standard test conditions (°C)
TEMP_COEFF_POWER = -0.004  # relative power loss per °C above STC for crystalline silicon (1/°C)
NOCT_TEMP_RISE = 0.03  # steady-state cell-temp rise per W/m² of POA (°C·m²/W)


class DifferentiableSolarPlant(nn.Module):
    """Differentiable physical model of a single metered solar site.

    Each physical parameter is a mean-field variational posterior N(mu, sigma^2).
    Training minimises `negative_elbo`: a Gaussian reconstruction likelihood (learnable
    observation noise) plus a KL term against a fixed prior (see `kl_divergence`). The KL
    term is what stops the posterior spreads from collapsing to zero under a pure
    reconstruction loss.

    Azimuth convention: measured from due south, east negative, west positive.
    """

    def __init__(
        self,
        prior_tilt: float,
        prior_azimuth: float,
        prior_dc_capacity: float,
        prior_ac_capacity: float,
    ) -> None:
        super().__init__()
        # Variational posteriors, each parameterised in an unconstrained space so no
        # gradient-breaking clamp is needed and the Normal posterior/prior stay
        # self-consistent: tilt in logit-space (squashed to (0, pi/2) in `forward`),
        # azimuth in radians, capacities in log-space (strictly positive).
        self.raw_tilt_mu = nn.Parameter(
            torch.special.logit(torch.tensor(prior_tilt) / (torch.pi / 2))
        )
        self.raw_tilt_log_std = nn.Parameter(torch.tensor(-2.0))

        self.azimuth_mu = nn.Parameter(torch.tensor(prior_azimuth))
        self.azimuth_log_std = nn.Parameter(torch.tensor(-2.0))

        self.log_dc_capacity_mu = nn.Parameter(torch.log(torch.tensor(prior_dc_capacity)))
        self.log_dc_capacity_log_std = nn.Parameter(torch.tensor(-2.0))

        self.log_ac_capacity_mu = nn.Parameter(torch.log(torch.tensor(prior_ac_capacity)))
        self.log_ac_capacity_log_std = nn.Parameter(torch.tensor(-2.0))

        # Observation-noise scale for the Gaussian reconstruction likelihood in
        # `negative_elbo`: learnable, in log-space to stay strictly positive.
        self.log_obs_noise = nn.Parameter(torch.tensor(0.0))

        # Fixed priors. Weakly-informative: "panels point roughly south at a UK roof pitch,
        # with a capacity near nameplate". These also keep the posterior spreads from collapsing.
        self.priors = {
            "raw_tilt": dist.Normal(
                torch.special.logit(torch.tensor(prior_tilt) / (torch.pi / 2)), torch.tensor(0.3)
            ),
            "azimuth": dist.Normal(torch.tensor(prior_azimuth), torch.tensor(0.5)),
            "log_dc_capacity": dist.Normal(
                torch.log(torch.tensor(prior_dc_capacity)), torch.tensor(0.5)
            ),
            "log_ac_capacity": dist.Normal(
                torch.log(torch.tensor(prior_ac_capacity)), torch.tensor(0.3)
            ),
        }

    def posteriors(self) -> dict[str, dist.Normal]:
        """The current variational posterior for each physical parameter."""
        return {
            "raw_tilt": dist.Normal(self.raw_tilt_mu, self.raw_tilt_log_std.exp()),
            "azimuth": dist.Normal(self.azimuth_mu, self.azimuth_log_std.exp()),
            "log_dc_capacity": dist.Normal(
                self.log_dc_capacity_mu, self.log_dc_capacity_log_std.exp()
            ),
            "log_ac_capacity": dist.Normal(
                self.log_ac_capacity_mu, self.log_ac_capacity_log_std.exp()
            ),
        }

    def kl_divergence(self) -> torch.Tensor:
        """Sum of KL(posterior || prior) over all parameters — the regulariser in the ELBO."""
        q, p = self.posteriors(), self.priors
        return sum(dist.kl_divergence(q[key], p[key]) for key in q)

    def negative_elbo(
        self, predicted_power: torch.Tensor, observed_power: torch.Tensor, dataset_size: int
    ) -> torch.Tensor:
        """Minibatch training loss: Gaussian reconstruction likelihood + scaled KL.

        The KL term regularises the whole-dataset objective once, so on a minibatch it is
        weighted by batch_size / dataset_size — otherwise the effective prior strength
        would depend on how the data happens to be batched.
        """
        obs_noise = self.log_obs_noise.exp()
        log_lik = dist.Normal(predicted_power, obs_noise).log_prob(observed_power).sum()
        kl_weight = observed_power.numel() / dataset_size
        return kl_weight * self.kl_divergence() - log_lik

    def forward(
        self,
        dni: torch.Tensor,  # direct normal irradiance (W/m^2)
        dhi: torch.Tensor,  # diffuse horizontal irradiance (W/m^2)
        ghi: torch.Tensor,  # global horizontal irradiance (W/m^2)
        sun_zenith_rad: torch.Tensor,
        sun_azimuth_rad: torch.Tensor,
        air_temp: torch.Tensor,  # ambient air temperature (°C)
    ) -> torch.Tensor:
        """Predict DC power with steady-state temperature derate. All inputs and operations preserve gradients."""
        q = self.posteriors()

        # Reparameterised samples (one per forward pass). tilt is sampled in logit-space
        # and squashed to (0, pi/2) — a sigmoid, not a clamp, so gradients flow everywhere
        # and the sample stays consistent with the Normal posterior used in the KL.
        tilt = torch.sigmoid(q["raw_tilt"].rsample()) * (torch.pi / 2)
        azimuth = q["azimuth"].rsample()
        dc_capacity = q["log_dc_capacity"].rsample().exp()  # strictly positive by construction

        # Angle of incidence on the tilted plane.
        cos_aoi = (
            torch.cos(sun_zenith_rad) * torch.cos(tilt)
            + torch.sin(sun_zenith_rad) * torch.sin(tilt) * torch.cos(sun_azimuth_rad - azimuth)
        ).clamp(min=0.0)  # no beam contribution when the sun is behind the panel

        # Isotropic-sky transposition to plane-of-array (POA) irradiance.
        poa_beam = dni * cos_aoi
        poa_sky_diffuse = dhi * (1.0 + torch.cos(tilt)) / 2.0
        poa_ground = ghi * GROUND_ALBEDO * (1.0 - torch.cos(tilt)) / 2.0
        poa = poa_beam + poa_sky_diffuse + poa_ground

        # DC power: capacity is the nameplate at the reference irradiance, derated for cell temperature.
        cell_temp = air_temp + NOCT_TEMP_RISE * poa
        temp_derate = 1.0 + TEMP_COEFF_POWER * (cell_temp - STC_CELL_TEMP)
        dc_power = dc_capacity * poa / STC_IRRADIANCE * temp_derate
        ac_capacity = q["log_ac_capacity"].rsample().exp()
        return torch.minimum(dc_power, ac_capacity)

Three physics details the sketch gets right:

  • Irradiance transposition. Plane-of-array (POA) irradiance is not GHI scaled by the angle of incidence — GHI already bakes in a cosine-of-zenith projection. We decompose the resource into beam, sky-diffuse, and ground-reflected components (an isotropic sky model) and transpose each correctly. The beam term uses DNI (not GHI) projected by the angle of incidence.
  • DC and AC capacity. The learnable dc_capacity is the DC nameplate in power units; POA is normalised by the reference irradiance (1000 W/m²) so that capacity falls out in MW at standard test conditions. ac_capacity is a separate learnable parameter that clips the inverter output via torch.minimum — a single plant clips hard, unlike the fleet soft clip.
  • Panel temperature derate. Cell temperature sits above ambient in proportion to absorbed POA irradiance; efficiency then falls roughly linearly with temperature above 25 °C. The sketch adds a steady-state derate using two module-level constants; in the full variational model these become learnable posteriors — and the subsection below extends this to the broken-cloud effect.

Angle convention: azimuth is measured from due south, with east negative and west positive (east = −90°, south = 0°, west = +90°). Beware that this differs from pvlib's convention (north = 0°, clockwise, in degrees); pvlib-pytorch should adopt pvlib's own convention and convert at the boundary — mismatched angle conventions are the classic silent bug in PV modelling.

Every fitted parameter is constrained by construction to a physically possible range, and that constraint does a different job from the prior. A tilt cannot leave 0° to 90°, a capacity cannot go negative, and an efficiency cannot exceed unity, because each posterior is parameterised in an unconstrained space and squashed into its admissible range. That range comes from a sigmoid or an exponential rather than a clamp, so gradients keep flowing everywhere. The constraint rules out what the physics forbids; the prior shapes what is merely implausible inside the range the constraint leaves. Prior widths are therefore a hyperparameter rather than a fixed choice. Where the experiment budget allows, we sweep them and read the answer off the forecast score.

The fitted parameters are effective parameters, not the plant's true parameters, and the evaluation follows from that. The energy-forecasting review reports Saint-Drenan et al. (2015) finding a fitted azimuth 5° off the surveyed value that simulated better than the surveyed value itself, because the fit absorbs the physical model's own systematic error. Saint-Drenan et al. concluded that the output "should be seen as a set of parameters that lead to the best simulation and not necessarily as the actual characteristics of the PV plant". Comparing a fitted tilt against a surveyed tilt is therefore a diagnostic — it catches a fit that has wandered somewhere physically absurd — and never the test. The test is the forecast score. The priors serve that same end: keeping the posterior spreads honest and regularising sites with little data, not pinning a posterior to a survey. That is why the priors in the sketch above are weakly informative rather than tight. No published work has run that test: the energy-forecasting review found no evidence for how much a photovoltaic power forecast improves from inferring tilt and azimuth rather than using a registered or nameplate value. The gain is therefore a hypothesis to test against the forecast score rather than a settled prize.

Fitting capacity itself as a distribution has a published precedent for wind, with a caveat on the headline number — the same treatment log_dc_capacity and log_ac_capacity get here as posteriors. The energy-forecasting review reports Pierrot and Pinson (2024)'s 34.2% CRPS gain from jointly fitting a time-varying capacity bound, and their own isolated test showing most of that gain comes from the joint fit rather than the bound alone.

Getting the model chain right is not a rounding error: the energy-forecasting review reports Mayer and Gróf (2021) finding a 13% mean-absolute-error gap between the best and worst of 32,400 model-chain combinations, even with every plant's geometry known exactly from its design documentation.

One detail the sketch deliberately ignores — interval averaging. NWP and satellite irradiance are half-hourly period means (period-ending in our NWP schema), but the sketch evaluates the solar geometry at an instant. Over 30 minutes the sun's hour angle moves ~7.5°, and near sunrise/sunset the transposition is strongly non-linear, so evaluating at a single timestamp biases the fit. Worse, a timestamp-convention error (period-start vs period-end vs mid-point) masquerades as an azimuth shift that the optimiser will happily absorb into the fitted parameters. The fix is cheap: evaluate the geometry at ~5-minute sub-steps across each half-hour and average the modelled power to the metered interval — and unit-test the timestamp convention end-to-end before trusting any fitted azimuth. This is not a V2-only concern: the same confusion is live in today's half-hourly NWP resample, where treating a period-ending radiation value as instantaneous shifts the modelled solar peak by up to three hours beyond forecast day 6. See NWP variable conventions, and the fix that lands first in v0.5.

Panel temperature

Steady-state derate. Cell temperature is not directly measured, and it is not a simple function of air temperature. On a hot, still, clear-sky summer day a panel can reach 60–70 °C — hot enough for efficiency to fall noticeably below its standard-test-condition (STC) value. What drives the cell above ambient is absorbed irradiance, moderated by convective cooling from wind. The standard Faiman relation captures this:

\[T_{\text{cell}} \approx T_{\text{air}} + \frac{\text{POA}}{U_0 + U_1\,v_{\text{wind}}}\]

Efficiency then falls roughly linearly with temperature above the 25 °C STC reference. The correction factor applied to DC output is:

\[\eta_T = 1 + \gamma\,(T_{\text{cell}} - 25\,^\circ\text{C})\]

with \(\gamma \approx -0.004\,/\,^\circ\text{C}\) (\(-0.4\,\%/^\circ\text{C}\)) for crystalline silicon. The constant NOCT_TEMP_RISE in the sketch is \(1/U_0\) with wind neglected — a useful simplification that removes the need for a wind-speed input in the code example. In the full variational model, \(U_0\), \(U_1\), and \(\gamma\) become learnable posteriors with tight physical priors, exactly like tilt, azimuth, and dc_capacity; the two new inputs (air temperature and POA) are already available: NWP temperature, and POA computed midway through the forward pass.

Thermal-mass upgrade and the broken-cloud effect. The steady-state model assumes the panel equilibrates to the current irradiance instantaneously. Real panels have thermal mass and lag by several minutes. This lag is exactly what produces the broken-cloud effect: a panel emerging cool from beneath cloud cover is briefly more efficient than one that has been baking under a clear sky, so the power peak immediately after cloud clearance can transiently exceed the steady-state clear-sky peak. This is also why peak daily yield on a partly-cloudy summer day can occasionally beat that on a fully clear day.

Capture this by promoting cell temperature from an instantaneous quantity to a dynamic latent state — first-order relaxation toward the equilibrium temperature with a learnable thermal time constant \(\tau\):

\[\frac{dT_{\text{cell}}}{dt} = \frac{T_{\text{eq}} - T_{\text{cell}}}{\tau}, \qquad T_{\text{eq}} = T_{\text{air}} + \frac{\text{POA}}{U_0 + U_1\,v_{\text{wind}}}\]

This is a state-space recurrence — the same pattern as a battery charge/discharge component.

Resolution caveat. The thermal time constant \(\tau\) is on the order of minutes, and the cloud transients that drive the broken-cloud effect occur on the same sub-half-hourly timescale. Half-hourly-mean POA discards exactly the sub-grid variability the dynamic model needs. This upgrade is therefore only worth making with higher-frequency irradiance inputs — satellite-derived irradiance at ~5-minute intervals, or sub-hourly metering — and the steady-state derate remains the right choice until that data is available.

Soiling

Soiling belongs in the v2 forward model, not in version 1. Version 1 reaches the same physics through a tabular stand-in: the drought and sustained-heat accumulator hands XGBoost the state variable \(d_t\) below — time since precipitation last exceeded a washing threshold — directly as a feature. This section makes the differentiable parameterisation concrete for when v2 needs it.

Capacity estimation names soiling as one of the shared physics corrections a learned fleet-scale model should absorb, and notes that a capacity estimate silently absorbs it otherwise. Soiling is worth modelling for Great Britain — not only for dustier climates. Dirt, dust, pollen, and bird droppings accumulate on panel glass and depress output by a few percent, and a decent fall of rain washes most of it off. Britain is normally rainy enough for the long-run average effect to be small, but the effect is driven by time since the last washing rainfall, not by climate averages, so a multi-month dry spell — London rooftops under Saharan dust after a rainless summer — is exactly the regime in which it stops being small. It is also one of the largest remaining unmodelled losses for any deployment in a dry or dusty climate.

How it would fit. Soiling is a multiplicative loss, so it enters as a ratio \(s_t \in (0, 1]\) applied to DC output alongside the temperature derate \(\eta_T\), driven by a small differentiable state that accumulates with dry time and resets with rain:

\[s_t = \max\!\left(s_{\min},\; 1 - \delta \cdot d_t\right), \qquad d_t = \begin{cases} 0 & \text{if } r_t > r_{\text{wash}} \\ d_{t-1} + \Delta t & \text{otherwise}\end{cases}\]

where \(d_t\) is time since the last washing rainfall and \(r_t\) is the rainfall rate. That is three learnable parameters — a soiling rate \(\delta\), a wash-off threshold \(r_{\text{wash}}\), and a floor \(s_{\min}\) — each with a tight physical prior, learned exactly like tilt and azimuth. The hard reset above wants softening (a sigmoid in \(r_t\), and a soft floor) so gradients flow, and the state-space recurrence is the same pattern as the thermal-mass upgrade above. No new input is needed: precipitation_surface is already in _ECMWF_ENS_VARS_TO_DOWNLOAD, so the rainfall history is in the NWP we already download.

Why it composes cleanly with the fleet node. UniversalSolarFleetNode carries magnitude in a non-decreasing per-week capacity series, because installed capacity only ever grows. Soiling is the opposite — a loss that reverses — so it cannot be expressed by bending that series. Trying to would destroy the monotone prior that is doing real work identifying installations. Factoring the two apart (monotone installed capacity × reversible soiling ratio) keeps both intact and adds only the three parameters above.

The real cost is identifiability, not implementation. Soiling and slower-than-assumed capacity growth both depress output, and they are separable only because soiling tracks rainfall history with a sawtooth shape whereas installations arrive as steps. Whether that separation actually holds at realistic noise levels needs demonstrating on synthetic data before the extra parameters are let anywhere near a real fit. In the aggregate-fleet case there is a second confound worth naming: real fleets are cleaned and rained on unevenly, so the fleet-level soiling ratio is a smeared average of many site-level sawtooths, which makes \(\delta\) easier to identify than \(r_{\text{wash}}\).

What a fitted five-parameter PV model measured on six solar farms

A throw-away experiment fitted this page's single-site physics on six metered solar farms without using any gradients, and three of its findings bear on the design sketched above. Panel tilt, panel azimuth, capacity, an inverter clipping limit, and a temperature coefficient were fitted per site by a Powell optimiser. The write-up is Does a weather product's beam/diffuse split help a PV forecast?, and its scope is six sites inside one 25 km by 23 km box in Lincolnshire over 2019 to 2026 — a measurement on this fleet rather than a general result about fitting PV models. The by-products that bear on estimating capacity rather than on this method, among them a half-hour timestamp error absorbed into the fitted azimuth, are on the capacity-estimation page.

A single site's five parameters are recoverable without gradients, so the case for differentiable physics cannot rest on that fit being hard. The objective is not convex. Refitting all 30 site-and-arm combinations from 64 independent random starting points each, a median of only 3 of the 64 starts reach the lowest loss found, and the worst start lands at up to 3.5 times that loss. What makes a derivative-free optimiser enough is the low dimension and a good fixed start rather than a single basin: the optimiser's own first start is a fixed vector shared by every seed, and that fixed start reaches the lowest loss in 19 of the 30 fits and is the only start to reach it in 6. The case for differentiable physics therefore rests on what this page argues elsewhere — posteriors that carry their own uncertainty, corrections shared across a whole fleet inside one training loop, gradients that reach back to the inputs, and continuity with the v2 engine — and not on a single site's PV parameters being hard to recover. None of that argues against differentiable physics; the finding narrows the list of reasons to reach for the method.

A five-parameter objective that is already multi-modal strengthens the case for a convex twin as the initialiser. The fleet node below carries far more parameters than five, and nothing in the single-site result suggests the fleet node's loss surface is better behaved. A starting point derived from a convex problem, rather than a hand-picked vector that happens to be good, is what that section proposes.

What limits the fitted physical model is its specification, not its optimiser. On the satellite irradiance source the fitted model reaches 6.28% of P99 output against XGBoost's 5.24%, so the tree is 1.05 percentage points better. A tree given nothing but the physical model's output lowers the error on neither irradiance source — the change is −0.025 percentage points [−0.061, +0.011] on the satellite source, an interval spanning zero. The deficit is therefore not a mis-calibration that rescaling the physical model's own output could repair.

The term that hurts is the transposition, and the sketch above carries no guard on it. The experiment compared a weather product's published direct-beam field against a beam estimated from the total irradiance by a separation model, and the fitted physical model and the tree disagree about which of the two beam fields is better. The physical model divides the horizontal beam by the cosine of the solar zenith angle to recover the direct normal irradiance, which magnifies a beam error without limit as the sun approaches the horizon: the physical model's preference between the two beam fields reaches +0.56 percentage points below 10 degrees of solar elevation, against +0.00 to +0.17 in the three elevation bands above it. Holding tilt and azimuth equal across the arms makes the disagreement larger rather than smaller, so the fitted geometry is not the cause — the write-up's comparison of the two instruments carries those numbers. DifferentiableSolarPlant takes the direct normal irradiance as an input, so the division sits upstream of the sketch wherever a weather product publishes beam on the horizontal plane, and the sketch's own clamp is on the cosine of the angle of incidence, which is a different quantity. The simplest guard is a floor on the zenith cosine, applied wherever the direct normal irradiance is derived.

A fitted parameter can land outside physics, and the fitted capacity then absorbs the slack. The fitted temperature coefficient runs from −0.0018 to +0.0020 per degree Celsius across the six sites, where a crystalline-silicon module's maximum-power coefficient is negative and datasheets cluster between −0.0045 and −0.0025 per degree Celsius — the range the sketch's fixed TEMP_COEFF_POWER of −0.004 sits inside. A positive fitted value means the parameter is absorbing something other than the temperature response, and the fitted capacity absorbs whatever derating the temperature coefficient leaves unapplied. The fitted tilts and azimuths are unaffected. A derivative-free fit with no prior has no way to rule a positive coefficient out, and the two devices the sketch above gives every parameter do: a parameterisation whose range excludes what the physics forbids, and a KL term pulling the posterior toward a physical prior. Making \(\gamma\) a learnable posterior, as Panel temperature proposes, is one place those two devices have concrete work to do.

Scaling to aggregate fleets: UniversalSolarFleetNode

A single tilt/azimuth pair cannot represent a "mishmash" of hundreds or thousands of rooftops. A fleet facing east, south, and west produces a broad, flat "mound" of power, whereas a single south-facing parameter produces a sharp "hill". Used on a primary substation, the single-array model will fail to fit the wide shoulders of morning and evening generation.

We fix this with a physics-informed basis expansion rather than simulating every individual system. (Aggregate, unmetered fleets behind a primary are a v2 concern — this node becomes one of the node types in the graph-structured engine.)

1. The basis-function insight

Think of the fleet not as a physical simulation of thousands of houses but as a signal-reconstruction problem. The aggregate curve of many fixed-tilt systems is a linear combination of a few fundamental shapes. East, south, and west span the space of fixed-tilt orientations: by mixing them (e.g. 30% east, 70% south) you can approximate almost any aggregate fixed-tilt curve, including SE/SW.

A basis mixture answers a documented failure of single-orientation fitting. The energy-forecasting review records that Saint-Drenan et al. (2015)'s algorithm "performs poorly" on an unmetered fleet, because it assumes one orientation per plant where the series is really the aggregated production of modules at many orientations. For a single metered site the same approach works well: Meng et al. (2020) infer the tilt and azimuth of 13 roof photovoltaic systems to mean absolute errors of 4.3° and 4.5° by matching the shape of metered output against irradiance measured elsewhere, needing no nameplate rating because both curves are normalised before matching. Meng et al. infer orientation and do not forecast power, though.

2. When to add a tracking basis

Single-axis trackers produce a flat-topped "top hat" shape that no mixture of fixed-tilt "bell curves" can reproduce. So where trackers are present, we add tracking as a fourth basis function. This matters mainly for commercial / utility ground-mount sitting behind a primary; domestic rooftop fleets are almost entirely fixed-tilt, so for a purely-domestic fleet the tracking weight should learn to ≈zero (or be dropped to save parameters).

3. Soft clipping for diverse inverters

A fleet contains many inverters of different sizes. An individual inverter hard-clips (a brick wall), but the aggregate clips softly: as irradiance rises the most undersized inverters saturate first, then the average inverters, then the oversized inverters — a smooth, curved shoulder rather than a sharp corner.

A tanh clip is tempting but wrong: P_max · tanh(P / P_max) already attenuates mid-range power (at P = P_max it returns only ≈0.76 P_max), so it biases the whole curve down, not just the shoulder. Instead use a smooth-min that stays roughly linear until the limit and only then rolls off:

\[\text{smin}(P, P_{\max}) = P_{\max} - \frac{1}{\beta}\,\operatorname{softplus}\!\bigl(\beta\,(P_{\max} - P)\bigr)\]

The learnable \(P_{\max}\) is the effective aggregate inverter (AC) capacity; the learnable sharpness \(\beta\) sets the curvature of the shoulder.

4. Implementation: UniversalSolarFleetNode

This upgrades the node to handle installed-capacity growth, orientation mix, tracking, and soft clipping in a single differentiable module. Note the division of labour: the softmax mix weights sum to 1, so the mixture alone is magnitude-free ("which way does the fleet face"); the magnitude ("how much is installed") is carried by a separate per-week capacity series, built as a cumulative sum of non-negative weekly increments so it is non-decreasing by construction — exactly the monotone representation from the disaggregation plan.

import torch
import torch.nn as nn
import torch.nn.functional as F


class UniversalSolarFleetNode(nn.Module):
    """Aggregate solar fleet behind one substation: capacity growth + orientation mix + soft clip."""

    def __init__(self, n_weeks: int) -> None:
        super().__init__()
        # Installed DC capacity, tracked per week as a cumulative sum of non-negative
        # increments (installs only ever add capacity — see the disaggregation
        # roadmap page). An L1 penalty on the increments (not shown) pushes most weeks
        # to exactly zero growth, because installs arrive in occasional bursts.
        self.raw_capacity_increments = nn.Parameter(torch.full((n_weeks,), -2.0))

        # Orientation mix over [east, south, west, tracking], tracked per week because the
        # fleet composition drifts over time as new systems are installed.
        self.mix_logits = nn.Parameter(torch.zeros(n_weeks, 4))

        # Effective aggregate inverter (AC) capacity and the curvature of the clip shoulder.
        # Both pushed through softplus so they stay strictly positive.
        self.raw_inverter_capacity = nn.Parameter(torch.tensor(2.3))  # softplus(2.3) ~ 10
        self.raw_clip_sharpness = nn.Parameter(torch.tensor(0.0))

    def forward(
        self, sun_vec: torch.Tensor, weather: torch.Tensor, week_idx: torch.Tensor
    ) -> torch.Tensor:
        # 1. The four physical basis curves (ideal shapes from geometry + weather), each
        #    normalised to unit DC capacity.
        #    calc_fixed / calc_tracker would come from pvlib-pytorch (see below).
        #    Azimuth convention matches the single-site model: south = 0, east = -90, west = +90.
        p_east = self.calc_fixed(sun_vec, weather, azimuth_deg=-90.0)
        p_south = self.calc_fixed(sun_vec, weather, azimuth_deg=0.0)
        p_west = self.calc_fixed(sun_vec, weather, azimuth_deg=90.0)
        p_track = self.calc_tracker(sun_vec, weather)  # single-axis backtracking logic

        # 2. Mix the unit bases with this week's weights (the "asset identification"),
        #    then scale by this week's installed capacity — non-decreasing by construction.
        basis = torch.stack([p_east, p_south, p_west, p_track], dim=-1)
        weights = torch.softmax(self.mix_logits[week_idx], dim=-1)
        capacity = F.softplus(self.raw_capacity_increments).cumsum(dim=0)
        p_raw = capacity[week_idx] * (basis * weights).sum(dim=-1)

        # 3. Aggregate soft clip. Smooth-min stays ~linear until p_raw nears the limit,
        #    then rolls off — unlike tanh it does not suppress mid-range power.
        limit = F.softplus(self.raw_inverter_capacity)
        beta = F.softplus(self.raw_clip_sharpness) + 1e-3
        return limit - F.softplus(beta * (limit - p_raw)) / beta

A convex twin, for initialisation and sanity checking

With the four basis curves held fixed, most of this node is fixed shapes × unknown coefficients: reparameterise (mix weights × capacity) as four non-negative per-basis capacity series, and the fleet's below-clip output is linear in them, with monotone growth and sparse increments available as hard convex constraints and an \(\ell_1\) penalty. The one part that breaks convexity is the soft clip, so the convex twin is exact only below the clip (or with the clip frozen). That is still enough for two jobs: a principled initialiser for the PyTorch fit, and an independent sanity check — if the convex twin and the trained node disagree materially on installed capacity, the neural fit deserves investigation before its answer is trusted.

pvlib-pytorch

pvlib-pytorch is a planned Open Climate Fix open-source library — a differentiable, PyTorch-native port of pvlib — that we intend to spin out of this project. It would generalise the hand-rolled transposition and panel geometry in the single-site and fleet sketches into a reusable, tested component; treat those sketches as the prototype it grows from.

Applications in this project