airpredict/packages/shared/src/prediction.ts

import { POLLUTANTS, Pollutant, StandardCode, STANDARD_LIMITS } from './pollutants.js';
import { PredictionRating } from './enums.js';

/** 单一污染物的散发模型参数 (per material, per pollutant) */
export interface EmissionParams {
  /** 最低平衡释放量 Y0 (mg/m²) */
  y0: number;
  /** 平衡释放量范围 Yp (mg/m²) */
  yp: number;
  /** 平衡释放量变化率 B (m²/m³) */
  b: number;
}

/** 一种材料在某空间的用量 + 各污染物散发参数 */
export interface SpaceMaterialInput {
  materialId: string;
  /** 使用面积 A (m²)。长度单位时调用方需先 ×0.1 折算 */
  usageAmount: number;
  params: Record<Pollutant, EmissionParams>;
}

/** 空间环境条件 */
export interface SpaceConditions {
  /** 体积 V (m³) */
  volume: number;
  /** 温度 T (°C) */
  temperature: number;
  /** 湿度 H (%rh，如 50；内部折算为小数) */
  humidity: number;
  /** 通风换气次数 ACH (次/小时) */
  ventilationRate: number;
  standard: StandardCode;
}

export interface MaterialContribution {
  materialId: string;
  /** 各污染物：该材料占最终浓度的份额 (mg/m³) = Gn × C */
  contribution: Record<Pollutant, number>;
  /** 各污染物贡献率 Gn (0~1) = (Yn修正·Ln)/C密闭 */
  contributionRate: Record<Pollutant, number>;
}

export interface SpacePredictionResult {
  /** 各污染物最终预测浓度 C (mg/m³) */
  concentration: Record<Pollutant, number>;
  /** 各污染物密闭浓度 C密闭 (mg/m³)，贡献率分母 */
  cClosed: Record<Pollutant, number>;
  /** 所用标准的各污染物限值 */
  limits: Record<Pollutant, number>;
  /** 各污染物是否超标 */
  exceeded: Record<Pollutant, boolean>;
  /** 各材料贡献 */
  contributions: MaterialContribution[];
  /** 各污染物的污染源材料ID（仅超标污染物；累计贡献率>50%的前几个材料） */
  sources: Record<Pollutant, string[]>;
  rating: PredictionRating;
}

/** 平衡释放量范围阈值：Yp < 0.01 视为该污染物不释放 */
const YP_THRESHOLD = 0.01;

const zero = (): Record<Pollutant, number> =>
  Object.fromEntries(POLLUTANTS.map((p) => [p, 0])) as Record<Pollutant, number>;

/**
 * 单空间多污染物预测（中国建材院模型 20240307）。
 * 每种污染物独立计算：
 *   Ln=A/V → 筛Yp≥0.01 → Ln修正=Σ(Li·Ypi)/Ypn → Yn修正=Y0+Yp·e^(−Ln修正/B)
 *   → C密闭=Σ(Yi修正·Li) → E=0.51+8.74·C密闭/Σ(Li·Ypi) → U=1/(1+E·ACH)
 *   → C通风=U·C密闭 → WS=e^(0.102(T−23))·e^(1.2312(H−0.5)) → C=WS·C通风
 *   贡献率 Gn=(Yn修正·Ln)/C密闭；污染源=超标污染物中按Gn降序累计>50%的材料。
 */
export function predictSpace(
  materials: SpaceMaterialInput[],
  cond: SpaceConditions,
): SpacePredictionResult {
  const V = Math.max(1e-6, cond.volume);
  const ACH = Math.max(0, cond.ventilationRate);
  const T = cond.temperature;
  const H = cond.humidity > 1 ? cond.humidity / 100 : cond.humidity; // %rh → 小数
  const WS = Math.exp(0.102 * (T - 23)) * Math.exp(1.2312 * (H - 0.5));
  const limits = STANDARD_LIMITS[cond.standard];

  // 真实承载率 Ln = A/V
  const Ln = new Map<string, number>();
  for (const m of materials) Ln.set(m.materialId, m.usageAmount / V);

  const concentration = zero();
  const cClosed = zero();
  const exceeded = {} as Record<Pollutant, boolean>;
  const sources = {} as Record<Pollutant, string[]>;
  const contribByMat: Record<string, Record<Pollutant, number>> = {};
  const rateByMat: Record<string, Record<Pollutant, number>> = {};
  for (const m of materials) {
    contribByMat[m.materialId] = zero();
    rateByMat[m.materialId] = zero();
  }

  for (const p of POLLUTANTS) {
    // 1. 筛选释放该污染物的材料（Yp≥0.01）
    const emitters = materials.filter((m) => (m.params[p]?.yp ?? 0) >= YP_THRESHOLD);
    // Σ(Li·Ypi)
    const sumLYp = emitters.reduce((s, m) => s + Ln.get(m.materialId)! * m.params[p].yp, 0);

    // 2. Yn修正 = Y0 + Yp·e^(−Ln修正/B)，其中 Ln修正 = Σ(Li·Ypi)/Ypn
    const yCorr = new Map<string, number>();
    for (const m of emitters) {
      const { y0, yp, b } = m.params[p];
      const lnCorr = yp > 0 ? sumLYp / yp : 0;
      const term = b > 0 ? yp * Math.exp(-lnCorr / b) : 0;
      yCorr.set(m.materialId, y0 + term);
    }

    // 3. 密闭浓度 C密闭 = Σ(Yn修正·Ln)
    const cClose = emitters.reduce((s, m) => s + yCorr.get(m.materialId)! * Ln.get(m.materialId)!, 0);
    cClosed[p] = round3(cClose);

    // 4. 通风修正
    const E = sumLYp > 0 ? 0.51 + 8.74 * (cClose / sumLYp) : 0;
    const U = 1 / (1 + E * ACH);
    const cVent = U * cClose;

    // 5. 温湿度修正 → 最终浓度
    const C = WS * cVent;
    concentration[p] = round3(C);

    // 6. 贡献率 Gn = (Yn修正·Ln)/C密闭
    for (const m of emitters) {
      const gn = cClose > 0 ? (yCorr.get(m.materialId)! * Ln.get(m.materialId)!) / cClose : 0;
      rateByMat[m.materialId][p] = gn;
      contribByMat[m.materialId][p] = round3(gn * C);
    }

    // 7. 超标判定 + 污染源（累计Gn>50%）
    const over = C > limits[p];
    exceeded[p] = over;
    if (over) {
      const sorted = emitters
        .map((m) => ({ id: m.materialId, gn: rateByMat[m.materialId][p] }))
        .sort((a, b) => b.gn - a.gn);
      const src: string[] = [];
      let cum = 0;
      for (const s of sorted) {
        src.push(s.id);
        cum += s.gn;
        if (cum > 0.5) break;
      }
      sources[p] = src;
    } else {
      sources[p] = [];
    }
  }

  const contributions: MaterialContribution[] = materials.map((m) => ({
    materialId: m.materialId,
    contribution: contribByMat[m.materialId],
    contributionRate: rateByMat[m.materialId],
  }));

  return {
    concentration,
    cClosed,
    limits,
    exceeded,
    contributions,
    sources,
    rating: ratingFor(concentration, cond.standard),
  };
}

function round3(n: number): number {
  return Math.round(n * 1000) / 1000;
}

/**
 * 评级（占位：官方评分/评级规则文档"待补充"）。
 * 现按各污染物相对限值的最差比值给级：A≤60% B≤100% C≤150% D>150%。
 */
export function ratingFor(
  concentration: Record<Pollutant, number>,
  standard: StandardCode,
): PredictionRating {
  const limits = STANDARD_LIMITS[standard];
  let worst = 0;
  for (const p of POLLUTANTS) {
    worst = Math.max(worst, limits[p] > 0 ? concentration[p] / limits[p] : 0);
  }
  if (worst <= 0.6) return 'A';
  if (worst <= 1.0) return 'B';
  if (worst <= 1.5) return 'C';
  return 'D';
}