Solar MPPT Power System

#include "plant.h"
#include <math.h>

/* Shape factor for the exponential I-V knee. 0.055 puts the maximum power
 * point at roughly 0.8*Voc, which is where a real silicon panel's sits. */
#define KNEE 0.055f

void plant_init(plant_t *p) {
    p->voc = 21.0f;             /* a 30W 12V-class panel in full sun */
    p->isc = 1.9f;
    p->batt_ocv = 12.0f;        /* a fairly discharged 12V SLA */
    p->batt_rint = 0.35f;
    p->batt_capacity_ah = 7.0f;
    p->efficiency = 0.92f;
    p->load_a = 0.0f;
    p->mppt_pgood = 1;
    p->bias_pgood = 1;
    p->vbus = p->batt_ocv;
    p->ichg = p->iin = p->ppanel = p->panel_v = 0.0f;
}

static float panel_current_at(const plant_t *p, float v) {
    if (v <= 0.0f) return p->isc;
    if (v >= p->voc) return 0.0f;
    float i = p->isc * (1.0f - expf((v - p->voc) / (KNEE * p->voc)));
    return i < 0.0f ? 0.0f : i;
}

float plant_panel_power_at(const plant_t *p, float v) {
    return v * panel_current_at(p, v);
}

float plant_panel_pmax(const plant_t *p, float *v_at_max) {
    float best = 0.0f, bv = 0.0f;
    for (int k = 1; k < 2000; k++) {
        float v = p->voc * (float)k / 2000.0f;
        float pw = plant_panel_power_at(p, v);
        if (pw > best) { best = pw; bv = v; }
    }
    if (v_at_max) *v_at_max = bv;
    return best;
}

/* Operating voltage on the HIGH-voltage branch of the I-V curve that
 * delivers exactly p_want, or the MPP voltage if p_want is unreachable.
 * The converter naturally settles on the high-voltage branch because that
 * is the stable side for a constant-power load. */
/* Overload collapse steepness: how fast delivered power falls once the
 * converter demands more than the panel's maximum. A real panel slides down
 * the low-voltage branch of its I-V curve, where power DECREASES as current
 * is drawn harder -- pinning it at Pmax instead would mean demanding too
 * much costs nothing, and MPPT would have no job to do. */
#define COLLAPSE_K 1.5f

static float panel_v_for_power(const plant_t *p, float p_want, float *p_got) {
    float vmp, pmax = plant_panel_pmax(p, &vmp);
    if (p_want >= pmax) {
        float overload = p_want / pmax - 1.0f;
        *p_got = pmax / (1.0f + COLLAPSE_K * overload);
        /* operating point has slid below Vmp, onto the low-voltage branch */
        return vmp * (*p_got / pmax);
    }

    float lo = vmp, hi = p->voc;           /* power decreases from vmp to voc */
    for (int k = 0; k < 60; k++) {
        float mid = 0.5f * (lo + hi);
        if (plant_panel_power_at(p, mid) > p_want) lo = mid; else hi = mid;
    }
    *p_got = plant_panel_power_at(p, 0.5f * (lo + hi));
    return 0.5f * (lo + hi);
}

void plant_step(plant_t *p, float vbatt_cmd, uint32_t dt_ms) {
    /* What the converter would deliver if the panel could support it: the
     * battery accepts (Vbus - OCV)/Rint, and the load takes its own share. */
    float ichg_want = (vbatt_cmd - p->batt_ocv) / p->batt_rint;
    if (ichg_want < 0.0f) ichg_want = 0.0f;          /* cannot sink current */
    float iout_want = ichg_want + p->load_a;
    float pout_want = vbatt_cmd * iout_want;
    float pin_want = pout_want / p->efficiency;

    float pin_got;
    p->panel_v = panel_v_for_power(p, pin_want, &pin_got);
    p->ppanel = pin_got;
    p->iin = (p->panel_v > 0.0f) ? pin_got / p->panel_v : 0.0f;

    float pout_got = pin_got * p->efficiency;

    if (pin_got >= pin_want - 1e-4f) {
        /* panel keeps up: the converter holds its commanded voltage */
        p->vbus = vbatt_cmd;
        p->ichg = ichg_want;
    } else {
        /* panel is the limit: the bus sags until supply and demand agree.
         * Solve Vbus*((Vbus-OCV)/Rint + load) = pout_got for Vbus. */
        float lo = p->batt_ocv, hi = vbatt_cmd;
        for (int k = 0; k < 60; k++) {
            float mid = 0.5f * (lo + hi);
            float pm = mid * ((mid - p->batt_ocv) / p->batt_rint + p->load_a);
            if (pm < pout_got) lo = mid; else hi = mid;
        }
        p->vbus = 0.5f * (lo + hi);
        p->ichg = (p->vbus - p->batt_ocv) / p->batt_rint;
        if (p->ichg < 0.0f) p->ichg = 0.0f;
    }

    /* the battery charges: OCV rises with accumulated charge */
    float hours = (float)dt_ms / 3600000.0f;
    float net_a = p->ichg - p->load_a;
    p->batt_ocv += (net_a * hours / p->batt_capacity_ah) * 2.4f;
    if (p->batt_ocv < 10.5f) p->batt_ocv = 10.5f;
    if (p->batt_ocv > 14.8f) p->batt_ocv = 14.8f;
}