machineGroupControl/src/optimizer/bepGravitation.js

// BEP-gravitation optimizer: bias flow allocation toward each pump's BEP,
// then refine via marginal-cost swaps. `ctx` shape matches bestCombination.js.

const MC_ITER_CAP = 50;        // marginal-cost refinement iterations
const MC_RELATIVE_EXIT = 0.001; // exit when the mc gap is < 0.1% of expensive.mc

// Estimate dP/dQ slopes around the BEP on the group operating point.
// Returns finite numbers for everything; falls back to zero slopes if the
// curve is flat or the machine has not been initialised.
function estimateSlopesAtBEP(machine, Q_BEP, ctx, delta = 1.0) {
    const { groupCurves } = ctx;
    const { groupFlow, groupNCog, groupCalcPower } = groupCurves;

    const minFlow = groupFlow(machine).currentFxyYMin;
    const maxFlow = groupFlow(machine).currentFxyYMax;
    const span = Math.max(0, maxFlow - minFlow);
    const normalizedCog = Math.max(0, Math.min(1, groupNCog(machine) || 0));
    const targetBEP = Q_BEP ?? (minFlow + span * normalizedCog);

    const clampFlow = (flow) => Math.min(maxFlow, Math.max(minFlow, flow));
    const center = clampFlow(targetBEP);
    const deltaSafe = Math.max(delta, 0.01);
    const leftFlow = clampFlow(center - deltaSafe);
    const rightFlow = clampFlow(center + deltaSafe);

    const powerAt = (flow) => groupCalcPower(machine, flow);
    const P_center = powerAt(center);
    const P_left = powerAt(leftFlow);
    const P_right = powerAt(rightFlow);
    const slopeLeft = (P_center - P_left) / Math.max(1e-6, center - leftFlow);
    const slopeRight = (P_right - P_center) / Math.max(1e-6, rightFlow - center);
    const alpha = Math.max(1e-6, (Math.abs(slopeLeft) + Math.abs(slopeRight)) / 2);

    return { slopeLeft, slopeRight, alpha, Q_BEP: center, P_BEP: P_center };
}

// Redistribute `delta` across pumps using slope-derived weights; flatter
// curves attract more flow. Bounded: exits on zero progress or no capacity.
function redistributeFlowBySlope(pumpInfos, flowDistribution, delta, directional = true) {
    const tolerance = 1e-3;
    let remaining = delta;
    const entryMap = new Map(flowDistribution.map(entry => [entry.machineId, entry]));

    while (Math.abs(remaining) > tolerance) {
        const increasing = remaining > 0;
        const candidates = pumpInfos.map(info => {
            const entry = entryMap.get(info.id);
            if (!entry) return null;
            const capacity = increasing ? info.maxFlow - entry.flow : entry.flow - info.minFlow;
            if (capacity <= tolerance) return null;
            const slope = increasing
                ? (directional ? info.slopes.slopeRight : info.slopes.alpha)
                : (directional ? info.slopes.slopeLeft : info.slopes.alpha);
            const weight = 1 / Math.max(1e-6, Math.abs(slope) || info.slopes.alpha || 1);
            return { entry, capacity, weight };
        }).filter(Boolean);

        if (!candidates.length) break;
        const weightSum = candidates.reduce((sum, c) => sum + c.weight * c.capacity, 0);
        if (weightSum <= 0) break;

        let progress = 0;
        candidates.forEach(candidate => {
            let share = (candidate.weight * candidate.capacity / weightSum) * Math.abs(remaining);
            share = Math.min(share, candidate.capacity);
            if (share <= 0) return;
            if (increasing) candidate.entry.flow += share;
            else candidate.entry.flow -= share;
            progress += share;
        });

        if (progress <= tolerance) break;
        remaining += increasing ? -progress : progress;
    }
}

function _marginalCostRefine(flowDistribution, pumpInfos, Qd, ctx) {
    const { groupCalcPower } = ctx.groupCurves;
    const mcDelta = Math.max(1e-6, (Qd / pumpInfos.length) * 0.005);

    for (let iter = 0; iter < MC_ITER_CAP; iter++) {
        const mcEntries = flowDistribution.map(entry => {
            const info = pumpInfos.find(i => i.id === entry.machineId);
            const pNow = groupCalcPower(info.machine, entry.flow);
            const pUp = groupCalcPower(info.machine, Math.min(info.maxFlow, entry.flow + mcDelta));
            return { entry, info, mc: (pUp - pNow) / mcDelta };
        });

        let expensive = null;
        let cheap = null;
        for (const e of mcEntries) {
            if (e.entry.flow > e.info.minFlow + mcDelta && (!expensive || e.mc > expensive.mc)) expensive = e;
            if (e.entry.flow < e.info.maxFlow - mcDelta && (!cheap || e.mc < cheap.mc)) cheap = e;
        }
        if (!expensive || !cheap || expensive === cheap) break;
        if (expensive.mc - cheap.mc < expensive.mc * MC_RELATIVE_EXIT) break;

        const before = groupCalcPower(expensive.info.machine, expensive.entry.flow)
                     + groupCalcPower(cheap.info.machine, cheap.entry.flow);
        const after = groupCalcPower(expensive.info.machine, expensive.entry.flow - mcDelta)
                    + groupCalcPower(cheap.info.machine, cheap.entry.flow + mcDelta);
        if (after < before) {
            expensive.entry.flow -= mcDelta;
            cheap.entry.flow += mcDelta;
        } else {
            break;
        }
    }
}

function calcBestCombinationBEPGravitation(combinations, Qd, ctx, method = 'BEP-Gravitation-Directional') {
    const { machines, groupCurves } = ctx;
    const { groupFlow, groupNCog, groupCalcPower } = groupCurves;
    const directional = method === 'BEP-Gravitation-Directional';

    let bestCombination = null;
    let bestPower = Infinity;
    let bestFlow = 0;
    let bestCog = 0;
    let bestDeviation = Infinity;

    combinations.forEach(combination => {
        const pumpInfos = combination.map(machineId => {
            const machine = machines[machineId];
            const minFlow = groupFlow(machine).currentFxyYMin;
            const maxFlow = groupFlow(machine).currentFxyYMax;
            const span = Math.max(0, maxFlow - minFlow);
            const NCog = Math.max(0, Math.min(1, groupNCog(machine) || 0));
            const estimatedBEP = minFlow + span * NCog;
            const slopes = estimateSlopesAtBEP(machine, estimatedBEP, ctx);
            return { id: machineId, machine, minFlow, maxFlow, NCog, Q_BEP: slopes.Q_BEP, slopes };
        });

        if (pumpInfos.length === 0) return;

        const flowDistribution = pumpInfos.map(info => ({
            machineId: info.id,
            flow: Math.min(info.maxFlow, Math.max(info.minFlow, info.Q_BEP)),
        }));

        let totalFlow = flowDistribution.reduce((s, e) => s + e.flow, 0);
        const delta = Qd - totalFlow;
        if (Math.abs(delta) > 1e-6) {
            redistributeFlowBySlope(pumpInfos, flowDistribution, delta, directional);
        }

        flowDistribution.forEach(entry => {
            const info = pumpInfos.find(i => i.id === entry.machineId);
            entry.flow = Math.min(info.maxFlow, Math.max(info.minFlow, entry.flow));
        });

        _marginalCostRefine(flowDistribution, pumpInfos, Qd, ctx);

        let totalPower = 0;
        totalFlow = 0;
        flowDistribution.forEach(entry => {
            totalFlow += entry.flow;
            const info = pumpInfos.find(i => i.id === entry.machineId);
            totalPower += groupCalcPower(info.machine, entry.flow);
        });

        const totalCog = pumpInfos.reduce((s, info) => s + info.NCog, 0);
        const deviation = pumpInfos.reduce((sum, info) => {
            const entry = flowDistribution.find(item => item.machineId === info.id);
            const deltaFlow = entry ? (entry.flow - info.Q_BEP) : 0;
            return sum + (deltaFlow * deltaFlow) * (info.slopes.alpha || 1);
        }, 0);

        const shouldUpdate = totalPower < bestPower
            || (totalPower === bestPower && deviation < bestDeviation);

        if (shouldUpdate) {
            bestCombination = flowDistribution.map(e => ({ ...e }));
            bestPower = totalPower;
            bestFlow = totalFlow;
            bestCog = totalCog;
            bestDeviation = deviation;
        }
    });

    return { bestCombination, bestPower, bestFlow, bestCog, bestDeviation, method };
}

module.exports = {
    calcBestCombinationBEPGravitation,
    estimateSlopesAtBEP,
    redistributeFlowBySlope,
};