machineGroupControl/test/integration/ncog-distribution.integration.test.js

/**
 * Group Distribution Strategy Comparison Test
 *
 * Compares three flow distribution strategies for a group of pumps:
 *   1. NCog/BEP-Gravitation (slope-weighted — favours pumps with flatter power curves)
 *   2. Equal distribution (same flow to every pump)
 *   3. Spillover (fill smallest pump first, overflow to next)
 *
 * For variable-speed centrifugal pumps, specific flow (Q/P) is monotonically
 * decreasing per pump (affinity laws: P ∝ Q³), so NCog = 0 for all pumps.
 * The real optimization value comes from the BEP-Gravitation algorithm's
 * slope-based redistribution, which IS sensitive to curve shape differences.
 *
 * These tests verify that:
 * - Asymmetric pumps produce different power slopes (the basis for optimization)
 * - BEP-Gravitation uses less total power than naive strategies for mixed pumps
 * - Equal pumps receive equal treatment under all strategies
 * - Spillover creates a visibly different distribution than BEP-weighted
 */
const test = require('node:test');
const assert = require('node:assert/strict');

const MachineGroup = require('../../src/specificClass');
const Machine = require('../../../rotatingMachine/src/specificClass');

const baseCurve = require('../../../generalFunctions/datasets/assetData/curves/hidrostal-H05K-S03R.json');

/* ---- helpers ---- */

function deepClone(obj) { return JSON.parse(JSON.stringify(obj)); }

function distortSeries(series, scale = 1, tilt = 0) {
  const last = series.length - 1;
  return series.map((v, i) => {
    const gradient = last === 0 ? 0 : i / last - 0.5;
    return Math.max(v * scale * (1 + tilt * gradient), 0);
  });
}

function createSyntheticCurve(mods) {
  const { flowScale = 1, powerScale = 1, flowTilt = 0, powerTilt = 0 } = mods;
  const curve = deepClone(baseCurve);
  Object.values(curve.nq).forEach(s => { s.y = distortSeries(s.y, flowScale, flowTilt); });
  Object.values(curve.np).forEach(s => { s.y = distortSeries(s.y, powerScale, powerTilt); });
  return curve;
}

const stateConfig = {
  time: { starting: 0, warmingup: 0, stopping: 0, coolingdown: 0 },
  movement: { speed: 1200, mode: 'staticspeed', maxSpeed: 1800 }
};

function createMachineConfig(id, label) {
  return {
    general: { logging: { enabled: false, logLevel: 'error' }, name: label, id, unit: 'm3/h' },
    functionality: { softwareType: 'machine', role: 'rotationaldevicecontroller' },
    asset: { model: 'hidrostal-H05K-S03R', unit: 'm3/h' },
    mode: {
      current: 'auto',
      allowedActions: { auto: ['execsequence', 'execmovement', 'flowmovement', 'statuscheck'] },
      allowedSources: { auto: ['parent', 'GUI'] }
    },
    sequences: {
      startup: ['starting', 'warmingup', 'operational'],
      shutdown: ['stopping', 'coolingdown', 'idle'],
      emergencystop: ['emergencystop', 'off'],
    }
  };
}

function createGroupConfig(name) {
  return {
    general: { logging: { enabled: false, logLevel: 'error' }, name },
    functionality: { softwareType: 'machinegroup', role: 'groupcontroller' },
    mode: { current: 'optimalcontrol' }
  };
}

/**
 * Bootstrap with differential pressure (upstream + downstream) so the predict
 * engine resolves a realistic fDimension and calcEfficiencyCurve produces
 * a proper BEP peak — not a monotonic Q/P curve.
 */
function bootstrapGroup(name, machineSpecs, diffMbar, upstreamMbar = 800) {
  const mg = new MachineGroup(createGroupConfig(name));
  const machines = {};
  for (const spec of machineSpecs) {
    const m = new Machine(createMachineConfig(spec.id, spec.label), stateConfig);
    if (spec.curveMods) m.updateCurve(createSyntheticCurve(spec.curveMods));
    // Set BOTH upstream and downstream so getMeasuredPressure computes differential
    m.updateMeasuredPressure(upstreamMbar, 'upstream', {
      timestamp: Date.now(), unit: 'mbar', childName: `pt-up-${spec.id}`, childId: `pt-up-${spec.id}`
    });
    m.updateMeasuredPressure(upstreamMbar + diffMbar, 'downstream', {
      timestamp: Date.now(), unit: 'mbar', childName: `pt-dn-${spec.id}`, childId: `pt-dn-${spec.id}`
    });
    mg.childRegistrationUtils.registerChild(m, 'downstream');
    machines[spec.id] = m;
  }
  return { mg, machines };
}

/** Distribute flow weighted by each machine's NCog (BEP position). */
function distributeByNCog(machines, Qd) {
  const entries = Object.entries(machines);
  let totalNCog = entries.reduce((s, [, m]) => s + (m.NCog || 0), 0);

  const distribution = {};
  for (const [id, m] of entries) {
    const min = m.predictFlow.currentFxyYMin;
    const max = m.predictFlow.currentFxyYMax;
    const flow = totalNCog > 0
      ? ((m.NCog || 0) / totalNCog) * Qd
      : Qd / entries.length;
    distribution[id] = Math.min(max, Math.max(min, flow));
  }

  let totalPower = 0;
  for (const [id, m] of entries) {
    totalPower += m.inputFlowCalcPower(distribution[id]);
  }
  return { distribution, totalPower };
}

/** Compute power at a given flow for a machine using its inverse curve. */
function powerAtFlow(machine, flow) {
  return machine.inputFlowCalcPower(flow);
}

/** Distribute by slope-weighting: flatter dP/dQ curves attract more flow. */
function distributeBySlopeWeight(machines, Qd) {
  const entries = Object.entries(machines);
  // Estimate slope (dP/dQ) at midpoint for each machine
  const pumpInfos = entries.map(([id, m]) => {
    const min = m.predictFlow.currentFxyYMin;
    const max = m.predictFlow.currentFxyYMax;
    const mid = (min + max) / 2;
    const delta = Math.max((max - min) * 0.05, 0.001);
    const pMid = powerAtFlow(m, mid);
    const pRight = powerAtFlow(m, Math.min(max, mid + delta));
    const slope = Math.abs((pRight - pMid) / delta);
    return { id, m, min, max, slope: Math.max(slope, 1e-6) };
  });

  // Weight = 1/slope: flatter curves get more flow
  const totalWeight = pumpInfos.reduce((s, p) => s + (1 / p.slope), 0);
  const distribution = {};
  let totalPower = 0;

  for (const p of pumpInfos) {
    const weight = (1 / p.slope) / totalWeight;
    const flow = Math.min(p.max, Math.max(p.min, Qd * weight));
    distribution[p.id] = flow;
    totalPower += powerAtFlow(p.m, flow);
  }

  return { distribution, totalPower };
}

/** Distribute equally. */
function distributeEqual(machines, Qd) {
  const entries = Object.entries(machines);
  const flowEach = Qd / entries.length;
  const distribution = {};
  let totalPower = 0;
  for (const [id, m] of entries) {
    const min = m.predictFlow.currentFxyYMin;
    const max = m.predictFlow.currentFxyYMax;
    const clamped = Math.min(max, Math.max(min, flowEach));
    distribution[id] = clamped;
    totalPower += powerAtFlow(m, clamped);
  }
  return { distribution, totalPower };
}

/** Spillover: fill smallest pump to max first, then overflow to next. */
function distributeSpillover(machines, Qd) {
  const entries = Object.entries(machines)
    .sort(([, a], [, b]) => a.predictFlow.currentFxyYMax - b.predictFlow.currentFxyYMax);
  let remaining = Qd;
  const distribution = {};
  let totalPower = 0;
  for (const [id, m] of entries) {
    const min = m.predictFlow.currentFxyYMin;
    const max = m.predictFlow.currentFxyYMax;
    const assigned = Math.min(max, Math.max(min, remaining));
    distribution[id] = assigned;
    remaining = Math.max(0, remaining - assigned);
  }
  for (const [id, m] of entries) {
    totalPower += powerAtFlow(m, distribution[id]);
  }
  return { distribution, totalPower };
}

/* ---- tests ---- */

test('NCog = 0 for centrifugal pumps (Q/P is monotonically decreasing with speed)', () => {
  // For variable-speed centrifugal pumps, P ∝ n³ and Q ∝ n, so Q/P ∝ 1/n²
  // which is always decreasing. Peak efficiency (Q/P) is always at index 0
  // (minimum speed), giving NCog = 0. This is physically correct — the MGC
  // compensates via slope-based redistribution instead.
  const { machines } = bootstrapGroup('ncog-basic', [
    { id: 'A', label: 'pump-A', curveMods: { flowScale: 1, powerScale: 1 } },
  ], 400); // 400 mbar differential

  const m = machines['A'];
  assert.ok(Number.isFinite(m.NCog), `NCog should be finite, got ${m.NCog}`);
  assert.strictEqual(m.NCog, 0, `NCog should be 0 for centrifugal pump (Q/P monotonically decreasing)`);
  assert.ok(m.cog > 0, `cog (peak specific flow) should be positive, got ${m.cog}`);
  assert.strictEqual(m.cogIndex, 0, `Peak Q/P should be at index 0 (minimum speed)`);
});

test('different curve shapes still yield NCog = 0 (Q/P limitation)', () => {
  // Even with powerTilt distortion, Q/P remains monotonically decreasing for
  // centrifugal pump curves because P grows much faster than Q with speed.
  // NCog = 0 for all shapes — the slope-based redistribution (tests 4-6)
  // is what actually differentiates asymmetric pumps.
  const { machines } = bootstrapGroup('ncog-shapes', [
    { id: 'early', label: 'early-BEP', curveMods: { flowScale: 1, powerScale: 1, powerTilt: 0.4 } },
    { id: 'late', label: 'late-BEP', curveMods: { flowScale: 1, powerScale: 1, powerTilt: -0.3 } },
  ], 400);

  const ncogEarly = machines['early'].NCog;
  const ncogLate = machines['late'].NCog;

  assert.strictEqual(ncogEarly, 0, `Early BEP NCog should be 0 (Q/P monotonic), got ${ncogEarly.toFixed(4)}`);
  assert.strictEqual(ncogLate, 0, `Late BEP NCog should be 0 (Q/P monotonic), got ${ncogLate.toFixed(4)}`);
  // Both cog values should still be computable and positive (peak Q/P at min speed)
  assert.ok(machines['early'].cog > 0, 'early cog should be positive');
  assert.ok(machines['late'].cog > 0, 'late cog should be positive');
});

test('NCog = 0 falls back to equal distribution (same as equal split)', () => {
  // When NCog = 0 for all pumps (centrifugal pump limitation), the
  // distributeByNCog helper falls back to equal distribution. This verifies
  // the fallback works correctly and produces the same result as explicit
  // equal distribution.
  const { machines } = bootstrapGroup('ncog-vs-equal', [
    { id: 'early', label: 'early-BEP', curveMods: { flowScale: 1, powerScale: 1, powerTilt: 0.4 } },
    { id: 'late', label: 'late-BEP', curveMods: { flowScale: 1, powerScale: 1, powerTilt: -0.3 } },
  ], 400);

  // Both NCog = 0 (confirmed by tests 1-2)
  assert.strictEqual(machines['early'].NCog, 0, 'early NCog should be 0');
  assert.strictEqual(machines['late'].NCog, 0, 'late NCog should be 0');

  const totalMax = machines['early'].predictFlow.currentFxyYMax + machines['late'].predictFlow.currentFxyYMax;
  const Qd = totalMax * 0.5;

  const ncogResult = distributeByNCog(machines, Qd);
  const equalResult = distributeEqual(machines, Qd);

  // With NCog = 0 for both, distributeByNCog falls back to equal split
  const ncogDiff = Math.abs(ncogResult.distribution['early'] - ncogResult.distribution['late']);
  const equalDiff = Math.abs(equalResult.distribution['early'] - equalResult.distribution['late']);
  assert.ok(
    Math.abs(ncogDiff - equalDiff) < Qd * 0.01,
    `NCog fallback should produce same distribution as equal split. ` +
    `ncogDiff=${ncogDiff.toFixed(4)}, equalDiff=${equalDiff.toFixed(4)}`
  );
});

test('asymmetric pumps have different power curve slopes', () => {
  // A pump with low powerScale has a flatter power curve
  const { machines } = bootstrapGroup('slope-check', [
    { id: 'flat', label: 'flat-power', curveMods: { flowScale: 1.2, powerScale: 0.7, flowTilt: 0.1 } },
    { id: 'steep', label: 'steep-power', curveMods: { flowScale: 0.8, powerScale: 1.4, flowTilt: -0.05 } },
  ], 400);

  // Compute slope at midpoint of each machine's range
  const slopes = {};
  for (const [id, m] of Object.entries(machines)) {
    const mid = (m.predictFlow.currentFxyYMin + m.predictFlow.currentFxyYMax) / 2;
    const delta = (m.predictFlow.currentFxyYMax - m.predictFlow.currentFxyYMin) * 0.05;
    const pMid = powerAtFlow(m, mid);
    const pRight = powerAtFlow(m, mid + delta);
    slopes[id] = (pRight - pMid) / delta;
  }

  assert.ok(slopes['flat'] > 0 && slopes['steep'] > 0, 'Both slopes should be positive');
  assert.ok(
    slopes['steep'] > slopes['flat'] * 1.3,
    `Steep pump should have notably higher slope. flat=${slopes['flat'].toFixed(0)}, steep=${slopes['steep'].toFixed(0)}`
  );
});

test('slope-weighted distribution routes more flow to flatter pump', () => {
  const { machines } = bootstrapGroup('slope-routing', [
    { id: 'flat', label: 'flat-power', curveMods: { flowScale: 1.2, powerScale: 0.7 } },
    { id: 'steep', label: 'steep-power', curveMods: { flowScale: 0.8, powerScale: 1.4 } },
  ], 400);

  const totalMax = machines['flat'].predictFlow.currentFxyYMax + machines['steep'].predictFlow.currentFxyYMax;
  const Qd = totalMax * 0.5;

  const slopeResult = distributeBySlopeWeight(machines, Qd);

  assert.ok(
    slopeResult.distribution['flat'] > slopeResult.distribution['steep'],
    `Flat pump should get more flow. flat=${slopeResult.distribution['flat'].toFixed(2)}, steep=${slopeResult.distribution['steep'].toFixed(2)}`
  );
});

test('slope-weighted uses less power than equal split for asymmetric pumps', () => {
  const { machines } = bootstrapGroup('power-compare', [
    { id: 'eff', label: 'efficient', curveMods: { flowScale: 1.2, powerScale: 0.7, flowTilt: 0.12 } },
    { id: 'std', label: 'standard', curveMods: { flowScale: 1, powerScale: 1 } },
  ], 400);

  const totalMax = machines['eff'].predictFlow.currentFxyYMax + machines['std'].predictFlow.currentFxyYMax;
  const demandLevels = [0.3, 0.5, 0.7].map(p => {
    const min = Math.max(machines['eff'].predictFlow.currentFxyYMin, machines['std'].predictFlow.currentFxyYMin);
    return min + (totalMax - min) * p;
  });

  let slopeWins = 0;
  const results = [];

  for (const Qd of demandLevels) {
    const slopeResult = distributeBySlopeWeight(machines, Qd);
    const equalResult = distributeEqual(machines, Qd);
    const spillResult = distributeSpillover(machines, Qd);

    results.push({
      demand: Qd,
      slopePower: slopeResult.totalPower,
      equalPower: equalResult.totalPower,
      spillPower: spillResult.totalPower,
    });

    if (slopeResult.totalPower <= equalResult.totalPower + 1) slopeWins++;
  }

  assert.ok(
    slopeWins >= 2,
    `Slope-weighted should use ≤ power than equal in ≥ 2/3 cases.\n` +
    results.map(r =>
      `  Qd=${r.demand.toFixed(1)}: slope=${r.slopePower.toFixed(1)}W, equal=${r.equalPower.toFixed(1)}W, spill=${r.spillPower.toFixed(1)}W`
    ).join('\n')
  );
});

test('spillover produces visibly different distribution than slope-weighted for mixed sizes', () => {
  const { machines } = bootstrapGroup('spillover-vs-slope', [
    { id: 'small', label: 'small-pump', curveMods: { flowScale: 0.6, powerScale: 0.55 } },
    { id: 'large', label: 'large-pump', curveMods: { flowScale: 1.5, powerScale: 1.2 } },
  ], 400);

  const totalMax = machines['small'].predictFlow.currentFxyYMax + machines['large'].predictFlow.currentFxyYMax;
  const Qd = totalMax * 0.5;

  const slopeResult = distributeBySlopeWeight(machines, Qd);
  const spillResult = distributeSpillover(machines, Qd);

  // Spillover fills the small pump first, slope-weight distributes by curve shape
  const slopeDiff = Math.abs(slopeResult.distribution['small'] - spillResult.distribution['small']);
  const percentDiff = (slopeDiff / Qd) * 100;

  assert.ok(
    percentDiff > 1,
    `Strategies should produce different distributions. ` +
    `Slope small=${slopeResult.distribution['small'].toFixed(2)}, ` +
    `Spill small=${spillResult.distribution['small'].toFixed(2)} (${percentDiff.toFixed(1)}% diff)`
  );
});

test('equal pumps get equal flow under all strategies', () => {
  const { machines } = bootstrapGroup('equal-pumps', [
    { id: 'A', label: 'pump-A', curveMods: { flowScale: 1, powerScale: 1 } },
    { id: 'B', label: 'pump-B', curveMods: { flowScale: 1, powerScale: 1 } },
  ], 400);

  const totalMax = machines['A'].predictFlow.currentFxyYMax + machines['B'].predictFlow.currentFxyYMax;
  const Qd = totalMax * 0.6;

  const slopeResult = distributeBySlopeWeight(machines, Qd);
  const equalResult = distributeEqual(machines, Qd);

  const tolerance = Qd * 0.01;

  assert.ok(
    Math.abs(slopeResult.distribution['A'] - slopeResult.distribution['B']) < tolerance,
    `Slope-weighted should split equally for identical pumps. A=${slopeResult.distribution['A'].toFixed(2)}, B=${slopeResult.distribution['B'].toFixed(2)}`
  );
  assert.ok(
    Math.abs(equalResult.distribution['A'] - equalResult.distribution['B']) < tolerance,
    `Equal should split equally. A=${equalResult.distribution['A'].toFixed(2)}, B=${equalResult.distribution['B'].toFixed(2)}`
  );

  // Power should be identical too
  assert.ok(
    Math.abs(slopeResult.totalPower - equalResult.totalPower) < 1,
    `Equal pumps should produce same total power under any strategy`
  );
});

test('full MGC optimalControl uses ≤ power than priorityControl for mixed pumps', async () => {
  const { mg, machines } = bootstrapGroup('mgc-full', [
    { id: 'eff', label: 'efficient', curveMods: { flowScale: 1.2, powerScale: 0.7, flowTilt: 0.1 } },
    { id: 'std', label: 'standard', curveMods: { flowScale: 1, powerScale: 1 } },
    { id: 'weak', label: 'weak', curveMods: { flowScale: 0.8, powerScale: 1.3, flowTilt: -0.08 } },
  ], 400);

  for (const m of Object.values(machines)) {
    await m.handleInput('parent', 'execSequence', 'startup');
  }

  // Run optimalControl. handleInput takes canonical m³/s post-refactor —
  // mirror the set.demand handler's percent → canonical mapping inline.
  mg.setMode('optimalcontrol');
  function pctCanonical(mgc, pct) {
    const dt = mgc.calcDynamicTotals();
    return mgc.interpolation.interpolate_lin_single_point(pct, 0, 100, dt.flow.min, dt.flow.max);
  }
  await mg.handleInput('parent', pctCanonical(mg, 50), Infinity);
  const optPower = mg.measurements.type('power').variant('predicted').position('atequipment').getCurrentValue() || 0;
  const optFlow = mg.measurements.type('flow').variant('predicted').position('atequipment').getCurrentValue() || 0;

  // Reset machines
  for (const m of Object.values(machines)) {
    await m.handleInput('parent', 'execSequence', 'shutdown');
    await m.handleInput('parent', 'execSequence', 'startup');
  }

  // Run priorityControl
  mg.setMode('prioritycontrol');
  await mg.handleInput('parent', pctCanonical(mg, 50), Infinity, ['eff', 'std', 'weak']);
  const prioPower = mg.measurements.type('power').variant('predicted').position('atequipment').getCurrentValue() || 0;
  const prioFlow = mg.measurements.type('flow').variant('predicted').position('atequipment').getCurrentValue() || 0;

  assert.ok(optFlow > 0, `Optimal should deliver flow, got ${optFlow}`);
  assert.ok(prioFlow > 0, `Priority should deliver flow, got ${prioFlow}`);

  // Compare efficiency (flow per unit power)
  const optEff = optPower > 0 ? optFlow / optPower : 0;
  const prioEff = prioPower > 0 ? prioFlow / prioPower : 0;

  assert.ok(
    optEff >= prioEff * 0.95,
    `Optimal efficiency should be ≥ priority (within 5% tolerance). ` +
    `Opt: ${optFlow.toFixed(1)}/${optPower.toFixed(1)}=${optEff.toFixed(6)} | ` +
    `Prio: ${prioFlow.toFixed(1)}/${prioPower.toFixed(1)}=${prioEff.toFixed(6)}`
  );
});