From 9b4b7b58b0b9358bbcb3b6e0eed655ce51adc454 Mon Sep 17 00:00:00 2001 From: sepehr Date: Sun, 19 Jul 2026 23:23:35 +0200 Subject: [PATCH] Stabilize Modelica Fixed/Free calibration and stop flaky Newton embeds. Sort constraint/control assembly for deterministic Jacobians, keep live z_ua on legacy HX, and add generic Fixed/Free UI assist without removing Z factors. Co-authored-by: Cursor --- .../src/components/panels/PropertiesPanel.tsx | 79 +++- apps/web/src/lib/calibAssist.test.ts | 226 ++++++++++ apps/web/src/lib/calibAssist.ts | 331 +++++++++++++++ apps/web/src/lib/configBuilder.test.ts | 47 +- apps/web/src/lib/configBuilder.ts | 51 ++- crates/cli/src/config.rs | 79 ++++ crates/cli/src/run.rs | 228 +++++++++- .../src/heat_exchanger/condenser.rs | 57 ++- .../src/heat_exchanger/evaporator.rs | 51 ++- crates/solver/src/inverse/bounded.rs | 2 +- crates/solver/src/inverse/constraint.rs | 2 +- crates/solver/src/inverse/embedding.rs | 28 +- crates/solver/src/inverse/expr.rs | 401 ++++++++++++++++++ crates/solver/src/inverse/mod.rs | 2 + crates/solver/src/system.rs | 226 +++++++++- 15 files changed, 1726 insertions(+), 84 deletions(-) create mode 100644 apps/web/src/lib/calibAssist.test.ts create mode 100644 apps/web/src/lib/calibAssist.ts create mode 100644 crates/solver/src/inverse/expr.rs diff --git a/apps/web/src/components/panels/PropertiesPanel.tsx b/apps/web/src/components/panels/PropertiesPanel.tsx index 1640bc5..80ce19c 100644 --- a/apps/web/src/components/panels/PropertiesPanel.tsx +++ b/apps/web/src/components/panels/PropertiesPanel.tsx @@ -39,6 +39,12 @@ import { } from "@/lib/configBuilder"; import { buildComponentInspector, getSolvedVariablesForComponent } from "@/lib/componentInspector"; import type { SolvedVariable } from "@/lib/api"; +import { + calibrateAdjacentActuatorPatches, + collectCalibHints, + isCalibRelevantNode, + withDerivedUaEff, +} from "@/lib/calibAssist"; import { ComponentIcon } from "@/components/canvas/ComponentIcon"; import ComponentDocPanel from "@/components/panels/ComponentDocPanel"; import { modelBannerForType } from "@/lib/componentDocMap"; @@ -138,11 +144,23 @@ export default function PropertiesPanel() { // Solver-computed unknowns owned by this component (free actuators + // calibration factors). Empty until a solve produces them. - const solvedVars = useMemo( - () => (node ? getSolvedVariablesForComponent(result, node.data.name) : []), - [node, result], + const solvedVars = useMemo(() => { + if (!node) return []; + const base = getSolvedVariablesForComponent(result, node.data.name); + return withDerivedUaEff(node.data.name, nodes, base); + }, [node, nodes, result]); + + const calibHints = useMemo( + () => (node ? collectCalibHints(nodes, edges, node) : []), + [node, nodes, edges], ); + /** One-click Free on adjacent actuators compatible with this Probe's Fixed measures. */ + const calibActuatorPatches = useMemo(() => { + if (!node || node.data.type !== "Probe") return []; + return calibrateAdjacentActuatorPatches(node, nodes, edges); + }, [node, nodes, edges]); + // After a successful solve, jump to Results when selecting a part (Dymola-like). useEffect(() => { if (!node) return; @@ -186,7 +204,8 @@ export default function PropertiesPanel() { Fixed ☐ = inconnue → il faut une équation de plus
  • - Calib : Probe Tsat Fixed (= équation) + Z_UA Free (= inconnue) + Calib : Probe mesure Fixed (= équation) + facteur{" "} + Free (= inconnue) — z_ua, z_dp, z_flow, f_w…
  • n_eq = n_unk. Pas de boucle de régulation.
  • @@ -305,6 +324,44 @@ export default function PropertiesPanel() { )} + {/* Always visible on Probe / Z-factor components (also in Results after solve). */} + {!isRegLoop && + isCalibRelevantNode(node) && + (calibHints.length > 0 || calibActuatorPatches.length > 0) && ( +
    +

    Fixed / Free (Modelica)

    + {calibHints.map((h, i) => ( +

    + {h.level === "ok" ? "✓ " : h.level === "warn" ? "⚠ " : "💡 "} + {h.message} +

    + ))} + {calibActuatorPatches.length > 0 && ( + + )} +
    + )} + {!isRegLoop && panelMode === "results" && ( {isRegLoop && (
    - Optionnel. Calibration SST + Z_UA → onglet Calibration du HX (cases - Fixed), pas ce nœud. + Optionnel. Calibration générique → onglet Calibration (Fixed/Free + sur mesures Probe + facteurs Z), pas ce nœud.
    )} @@ -548,12 +605,12 @@ export default function PropertiesPanel() { {activeTab === "Calibration" && (

    - Modelica : Probe Fixed → équation ( - Tsat = …). Facteur Free → inconnue ( - parameter z_ua(fixed=false)). Le start est la valeur affichée. + Fixed = parameter connu (équation si Probe).{" "} + Free = inconnue Newton. Pairing : Tsat/Tsh↔z_ua, + P↔z_dp, Capacity↔z_flow, T↔f_w. Tous les Z restent disponibles.

    -

    - Laisse UA override vide. n_eq doit égaler n_unk. +

    + Ex. UA nominale × Z_UA → UA_eff. Start = valeur affichée.

    {solvedVars.length > 0 && ( diff --git a/apps/web/src/lib/calibAssist.test.ts b/apps/web/src/lib/calibAssist.test.ts new file mode 100644 index 0000000..36e5374 --- /dev/null +++ b/apps/web/src/lib/calibAssist.test.ts @@ -0,0 +1,226 @@ +import { describe, it, expect } from "vitest"; +import type { Edge, Node } from "@xyflow/react"; +import { fixedFlagKey } from "./componentMeta"; +import { + adjacentProbeNames, + calibrateAdjacentActuatorPatches, + collectCalibHints, + withDerivedUaEff, + type CalibNodeData, +} from "./calibAssist"; +import type { SolvedVariable } from "./api"; + +function node( + id: string, + type: string, + name: string, + params: Record = {}, +): Node { + return { + id, + type: "entropykNode", + position: { x: 0, y: 0 }, + data: { type, name, circuit: 0, params }, + }; +} + +describe("calibAssist — Fixed/Free générique", () => { + it("lists adjacent Probe names for any component", () => { + const nodes = [ + node("e", "Evaporator", "evap", { ua: 6000, z_ua: 1 }), + node("p1", "Probe", "sst", { tsat_c: 5, [fixedFlagKey("tsat_c")]: true }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "e", target: "p1", sourceHandle: "outlet", targetHandle: "inlet" }, + ]; + expect(adjacentProbeNames("evap", nodes, edges)).toEqual(["sst"]); + }); + + it("one-click Free Z_UA when Probe Tsat Fixed", () => { + const nodes = [ + node("e", "FloodedEvaporator", "evap", { + ua: 9000, + z_ua: 1.0, + [fixedFlagKey("z_ua")]: true, + }), + node("p1", "Probe", "sst", { + tsat_c: 5, + [fixedFlagKey("tsat_c")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "e", target: "p1", sourceHandle: "outlet", targetHandle: "inlet" }, + ]; + const patches = calibrateAdjacentActuatorPatches(nodes[1], nodes, edges); + expect(patches).toHaveLength(1); + expect(patches[0]).toMatchObject({ + nodeId: "e", + params: { [fixedFlagKey("z_ua")]: false, z_ua: 1.0 }, + }); + }); + + it("one-click Free z_dp when Probe P Fixed", () => { + const nodes = [ + node("c", "Condenser", "cond", { + ua: 5000, + z_dp: 1.0, + [fixedFlagKey("z_dp")]: true, + z_ua: 1.0, + [fixedFlagKey("z_ua")]: true, + }), + node("p1", "Probe", "p_probe", { + p_bar: 12, + [fixedFlagKey("p_bar")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "c", target: "p1", sourceHandle: "inlet", targetHandle: "outlet" }, + ]; + const patches = calibrateAdjacentActuatorPatches(nodes[1], nodes, edges); + expect(patches).toHaveLength(1); + expect(patches[0].params[fixedFlagKey("z_dp")]).toBe(false); + // z_ua not compatible with P — stays Fixed + expect(patches[0].params[fixedFlagKey("z_ua")]).toBe(true); + }); + + it("one-click Free z_flow when Probe Capacity Fixed (compressor)", () => { + const nodes = [ + node("comp", "IsentropicCompressor", "comp", { + z_flow: 1.0, + [fixedFlagKey("z_flow")]: true, + f_w: 0.1, + [fixedFlagKey("f_w")]: true, + }), + node("p1", "Probe", "cap", { + capacity_w: 12000, + [fixedFlagKey("capacity_w")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "comp", target: "p1", sourceHandle: "outlet", targetHandle: "inlet" }, + ]; + const patches = calibrateAdjacentActuatorPatches(nodes[1], nodes, edges); + expect(patches).toHaveLength(1); + expect(patches[0].params[fixedFlagKey("z_flow")]).toBe(false); + expect(patches[0].params[fixedFlagKey("f_w")]).toBe(true); + }); + + it("one-click Free f_w when Probe T Fixed", () => { + const nodes = [ + node("comp", "IsentropicCompressor", "comp", { + f_w: 0.05, + [fixedFlagKey("f_w")]: true, + z_flow: 1.0, + [fixedFlagKey("z_flow")]: true, + }), + node("p1", "Probe", "dgt", { + t_c: 75, + [fixedFlagKey("t_c")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "comp", target: "p1", sourceHandle: "outlet", targetHandle: "inlet" }, + ]; + const patches = calibrateAdjacentActuatorPatches(nodes[1], nodes, edges); + expect(patches[0].params[fixedFlagKey("f_w")]).toBe(false); + expect(patches[0].params[fixedFlagKey("z_flow")]).toBe(true); + }); + + it("returns empty when compatible actuators already Free", () => { + const nodes = [ + node("e", "Evaporator", "evap", { + ua: 6000, + z_ua: 1.0, + [fixedFlagKey("z_ua")]: false, + }), + node("p1", "Probe", "sst", { + tsat_c: 5, + [fixedFlagKey("tsat_c")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "e", target: "p1", sourceHandle: "outlet", targetHandle: "inlet" }, + ]; + expect(calibrateAdjacentActuatorPatches(nodes[1], nodes, edges)).toHaveLength(0); + }); + + it("tips when Probe Fixed but compatible factor still Fixed", () => { + const nodes = [ + node("e", "Condenser", "cond", { + ua: 5000, + z_ua: 1.0, + [fixedFlagKey("z_ua")]: true, + }), + node("p1", "Probe", "sdt", { + tsat_c: 40, + [fixedFlagKey("tsat_c")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "e", target: "p1", sourceHandle: "inlet", targetHandle: "outlet" }, + ]; + const hints = collectCalibHints(nodes, edges, nodes[1]); + expect(hints.some((h) => h.level === "tip" && h.message.includes("décoche Fixed"))).toBe( + true, + ); + }); + + it("ok when Probe Fixed + compatible factor Free", () => { + const nodes = [ + node("e", "Condenser", "cond", { + ua: 5000, + z_dp: 1.0, + [fixedFlagKey("z_dp")]: false, + }), + node("p1", "Probe", "p_probe", { + p_bar: 12, + [fixedFlagKey("p_bar")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "e", target: "p1", sourceHandle: "inlet", targetHandle: "outlet" }, + ]; + const hints = collectCalibHints(nodes, edges, nodes[0]); + expect(hints.some((h) => h.level === "ok" && h.message.includes("z_dp"))).toBe(true); + }); + + it("warns Free factor without adjacent compatible Probe", () => { + const nodes = [ + node("e", "Evaporator", "evap", { + ua: 6000, + z_ua: 1.0, + [fixedFlagKey("z_ua")]: false, + }), + ]; + const hints = collectCalibHints(nodes, [], nodes[0]); + expect(hints.some((h) => h.level === "warn")).toBe(true); + }); + + it("always shows a baseline tip on Probe / HX when nothing paired yet", () => { + const probe = node("p1", "Probe", "sst", { tsat_c: 5 }); + const hx = node("e", "Evaporator", "evap", { ua: 6000, z_ua: 1.0 }); + expect(collectCalibHints([probe], [], probe).some((h) => h.level === "tip")).toBe(true); + expect(collectCalibHints([hx], [], hx).some((h) => h.level === "tip")).toBe(true); + }); + + it("appends derived UA_eff = UA × z_ua", () => { + const nodes = [node("e", "Evaporator", "evap", { ua: 8000, z_ua: 1 })]; + const solved: SolvedVariable[] = [ + { + id: "evap__z_ua", + component: "evap", + variable: "z_ua", + value: 0.85, + min: 0.05, + max: 3, + }, + ]; + const out = withDerivedUaEff("evap", nodes, solved); + expect(out).toHaveLength(2); + expect(out[1]).toMatchObject({ + variable: "UA_eff", + value: 6800, + component: "evap", + }); + }); +}); diff --git a/apps/web/src/lib/calibAssist.ts b/apps/web/src/lib/calibAssist.ts new file mode 100644 index 0000000..d262890 --- /dev/null +++ b/apps/web/src/lib/calibAssist.ts @@ -0,0 +1,331 @@ +/** + * Generic Modelica Fixed / Free (Unfixed) calibration assistants. + * + * Engine already emits embeddings for ANY Free `actuatorFactor` + Fixed Probe. + * This module only explains DoF balance and offers one-click Free on compatible + * adjacent actuators. Nothing is removed from the catalogue (all Z stay). + */ + +import type { Edge, Node } from "@xyflow/react"; +import { + COMPONENT_BY_TYPE, + fixedFlagKey, + isParamFixed, + type ParamMeta, +} from "./componentMeta"; +import { FACTOR_COMPATIBILITY, PROBE_PARAM_KIND } from "./configBuilder"; +import type { SolvedVariable } from "./api"; + +/** Minimal node data shape (avoids heavy coupling beyond pairing tables). */ +export type CalibNodeData = { + type: string; + name: string; + params: Record; + circuit?: number; +}; + +export type CalibHintLevel = "ok" | "tip" | "warn"; + +export interface CalibHint { + level: CalibHintLevel; + /** Short UI message (French). */ + message: string; + component?: string; +} + +export interface CalibParamPatch { + nodeId: string; + params: Record; +} + +type FixedProbeMeasure = { + key: string; + label: string; + kind: string; + factors: readonly string[]; +}; + +function actuatorParams(type: string): ParamMeta[] { + const meta = COMPONENT_BY_TYPE[type]; + if (!meta) return []; + return meta.params.filter((p) => p.fixable && p.actuatorFactor); +} + +function fixedProbeMeasures(probe: Node): FixedProbeMeasure[] { + const meta = COMPONENT_BY_TYPE.Probe; + if (!meta || probe.data.type !== "Probe") return []; + const params = probe.data.params ?? {}; + const out: FixedProbeMeasure[] = []; + + for (const p of meta.params) { + if (!p.fixable || !p.measureOutput) continue; + if (!isParamFixed(params, p)) continue; + const raw = params[p.key]; + if (raw === undefined || raw === "") continue; + const n = typeof raw === "number" ? raw : Number(raw); + if (!Number.isFinite(n)) continue; + const kind = PROBE_PARAM_KIND[p.key] ?? p.key; + const factors = FACTOR_COMPATIBILITY[kind] ?? []; + if (factors.length === 0) continue; + out.push({ key: p.key, label: p.label, kind, factors }); + } + + // Legacy measure+target + if (out.length === 0) { + const fixedFlag = params.__fixed_target; + if (fixedFlag === true || fixedFlag === "true") { + const kind = String(params.measure ?? ""); + const factors = FACTOR_COMPATIBILITY[kind] ?? []; + if (kind && params.target != null && factors.length > 0) { + out.push({ + key: "target", + label: kind, + kind, + factors, + }); + } + } + } + + return out; +} + +/** Neighbours connected by an edge (any component type). */ +export function adjacentNodes( + node: Node, + nodes: Node[], + edges: Edge[], +): Node[] { + const byId = new Map(nodes.map((n) => [n.id, n])); + const out: Node[] = []; + const seen = new Set(); + for (const e of edges) { + const otherId = e.source === node.id ? e.target : e.target === node.id ? e.source : null; + if (!otherId || seen.has(otherId)) continue; + const other = byId.get(otherId); + if (!other) continue; + seen.add(otherId); + out.push(other); + } + return out; +} + +/** Probe names adjacent to a named component. */ +export function adjacentProbeNames( + componentName: string, + nodes: Node[], + edges: Edge[], +): string[] { + const byName = new Map(nodes.map((n) => [n.data.name, n])); + const comp = byName.get(componentName); + if (!comp) return []; + return adjacentNodes(comp, nodes, edges) + .filter((n) => n.data.type === "Probe") + .map((n) => n.data.name); +} + +/** True if the node participates in Fixed/Free calibration UX. */ +export function isCalibRelevantNode(node: Node | null | undefined): boolean { + if (!node) return false; + if (node.data.type === "Probe") return true; + return actuatorParams(node.data.type).length > 0; +} + +/** + * Collect Modelica Fixed/Free tips for the selected node (or whole diagram). + * Does not replace hard Free-without-equation errors from validateConfig. + * Always returns at least a baseline tip for Probe / actuator components. + */ +export function collectCalibHints( + nodes: Node[], + edges: Edge[], + focus?: Node | null, +): CalibHint[] { + const hints: CalibHint[] = []; + + for (const probe of nodes) { + if (probe.data.type !== "Probe") continue; + const measures = fixedProbeMeasures(probe); + if (measures.length === 0) continue; + + const neighbours = adjacentNodes(probe, nodes, edges).filter( + (n) => n.data.type !== "Probe", + ); + + for (const measure of measures) { + for (const neigh of neighbours) { + if ( + focus && + focus.id !== probe.id && + focus.id !== neigh.id + ) { + continue; + } + for (const ap of actuatorParams(neigh.data.type)) { + const factor = ap.actuatorFactor!; + if (!measure.factors.includes(factor)) continue; + const free = !isParamFixed(neigh.data.params ?? {}, ap); + if (free) { + hints.push({ + level: "ok", + component: neigh.data.name, + message: + `OK Modelica : Probe « ${probe.data.name} ».${measure.label} Fixed ` + + `(équation) + ${factor} Free sur « ${neigh.data.name} » (inconnue).`, + }); + } else { + hints.push({ + level: "tip", + component: neigh.data.name, + message: + `Pour calibrer : décoche Fixed sur ${factor} de « ${neigh.data.name} » ` + + `(Probe « ${probe.data.name} ».${measure.label} impose déjà la mesure).`, + }); + } + } + } + } + } + + if (focus) { + for (const ap of actuatorParams(focus.data.type)) { + if (isParamFixed(focus.data.params ?? {}, ap)) continue; + const factor = ap.actuatorFactor!; + const probes = adjacentNodes(focus, nodes, edges).filter((n) => n.data.type === "Probe"); + const hasEquation = probes.some((pr) => + fixedProbeMeasures(pr).some((m) => m.factors.includes(factor)), + ); + if (!hasEquation) { + hints.push({ + level: "warn", + component: focus.data.name, + message: + `${factor} est Free (inconnue) : ajoute une Probe Fixed compatible ` + + `(${compatibleMeasureLabels(factor).join(", ") || "mesure"}) ` + + `sur le circuit, ou remets ${factor} Fixed.`, + }); + } + } + + // Baseline always visible — otherwise users see nothing when all Fixed (default). + if (hints.length === 0 && isCalibRelevantNode(focus)) { + if (focus.data.type === "Probe") { + const measures = fixedProbeMeasures(focus); + if (measures.length === 0) { + hints.push({ + level: "tip", + component: focus.data.name, + message: + "Fixed/Free : coche Fixed sur une mesure (Tsat, P, T, Capacity…) " + + "puis libère le facteur compatible sur le composant adjacent " + + "(z_ua, z_dp, z_flow, f_w).", + }); + } + } else { + const factors = actuatorParams(focus.data.type) + .map((p) => p.actuatorFactor) + .filter(Boolean) + .join(", "); + hints.push({ + level: "tip", + component: focus.data.name, + message: + `Fixed/Free : facteurs ${factors || "Z"} — Fixed = connu, Free = inconnue. ` + + `Pour calibrer : Probe Fixed adjacent + décoche Fixed sur le facteur.`, + }); + } + } + } + + return hints; +} + +function compatibleMeasureLabels(factor: string): string[] { + const labels: string[] = []; + for (const [kind, factors] of Object.entries(FACTOR_COMPATIBILITY)) { + if (factors.includes(factor)) labels.push(kind); + } + return labels; +} + +/** + * One-click: Free every Fixed actuator on neighbours that is compatible with + * this Probe's Fixed measures. Returns patches for `updateNodesParams`. + */ +export function calibrateAdjacentActuatorPatches( + probeNode: Node, + nodes: Node[], + edges: Edge[], +): CalibParamPatch[] { + if (probeNode.data.type !== "Probe") return []; + const measures = fixedProbeMeasures(probeNode); + if (measures.length === 0) return []; + + const wanted = new Set(measures.flatMap((m) => [...m.factors])); + const byId = new Map(); + + for (const neigh of adjacentNodes(probeNode, nodes, edges)) { + if (neigh.data.type === "Probe") continue; + let next = { ...(neigh.data.params ?? {}) }; + let changed = false; + for (const ap of actuatorParams(neigh.data.type)) { + const factor = ap.actuatorFactor!; + if (!wanted.has(factor)) continue; + if (!isParamFixed(next, ap)) continue; + next = { + ...next, + [fixedFlagKey(ap.key)]: false, + [ap.key]: + typeof next[ap.key] === "number" + ? next[ap.key] + : (ap.default as number | undefined) ?? 1.0, + }; + changed = true; + } + if (changed) { + byId.set(neigh.id, { nodeId: neigh.id, params: next }); + } + } + + return [...byId.values()]; +} + +/** @deprecated Prefer {@link calibrateAdjacentActuatorPatches}. */ +export function calibrateAdjacentZUaPatch( + probeNode: Node, + nodes: Node[], + edges: Edge[], +): CalibParamPatch | null { + const patches = calibrateAdjacentActuatorPatches(probeNode, nodes, edges); + const zUa = patches.find((p) => p.params[fixedFlagKey("z_ua")] === false); + return zUa ?? patches[0] ?? null; +} + +/** + * Append derived UA_eff = UA × z_ua when both are known after solve. + */ +export function withDerivedUaEff( + componentName: string, + nodes: Node[], + solved: SolvedVariable[], +): SolvedVariable[] { + const node = nodes.find((n) => n.data.name === componentName); + if (!node) return solved; + const z = solved.find((s) => s.variable === "z_ua"); + if (!z) return solved; + const uaRaw = node.data.params?.ua; + const ua = typeof uaRaw === "number" ? uaRaw : Number(uaRaw); + if (!Number.isFinite(ua) || ua <= 0) return solved; + if (solved.some((s) => s.variable === "UA_eff")) return solved; + return [ + ...solved, + { + id: `${componentName}__UA_eff`, + component: componentName, + variable: "UA_eff", + value: ua * z.value, + min: 0, + max: ua * z.max, + }, + ]; +} diff --git a/apps/web/src/lib/configBuilder.test.ts b/apps/web/src/lib/configBuilder.test.ts index ee70b24..6571d7d 100644 --- a/apps/web/src/lib/configBuilder.test.ts +++ b/apps/web/src/lib/configBuilder.test.ts @@ -505,7 +505,7 @@ describe("Modelica embeddings (Fixed / Free)", () => { expect(embeddings).toHaveLength(1); expect(embeddings[0]).toMatchObject({ id: "emb_evap_z_ua", - unknown: { component: "evap", factor: "z_ua", start: 0.3, min: 0.05 }, + unknown: { component: "evap", factor: "z_ua", start: 1.0, min: 0.05 }, equation: { component: "sst_probe", output: "saturationTemperature", @@ -519,7 +519,50 @@ describe("Modelica embeddings (Fixed / Free)", () => { const comp = cfg.circuits[0].components[0]; expect(comp[fixedFlagKey("z_ua")]).toBeUndefined(); expect(comp.z_ua).toBe(1.0); - expect(comp.ua).toBeUndefined(); // Free z_ua strips absolute ua + // Legacy HX: ua is the nominal base — keep it when Free z_ua. + expect(comp.ua).toBe(9000); + }); + + it("strips BPHX ua override when Z_UA is Free (geometry = UA_nominal)", () => { + const nodes = [ + node("c", "BphxCondenser", "cond", 0, { + n_plates: 40, + ua: 2500, + z_ua: 1.0, + [fixedFlagKey("z_ua")]: false, + }), + node("p1", "Probe", "sdt_probe", 0, { + fluid: "R134a", + tsat_c: 42.0, + [fixedFlagKey("tsat_c")]: true, + }), + ]; + const edges: Edge[] = [ + { id: "pe", source: "c", target: "p1", sourceHandle: "inlet", targetHandle: "outlet" }, + ]; + const cfg = buildScenarioConfig(nodes, edges); + expect(cfg.embeddings).toHaveLength(1); + expect(cfg.circuits[0].components[0].ua).toBeUndefined(); + }); + + it("fills missing required ua on FloodedEvaporator when Z_UA is Free", () => { + const nodes = [ + node("e", "FloodedEvaporator", "evap", 0, { + // ua intentionally missing / cleared in the canvas + z_ua: 1.0, + [fixedFlagKey("z_ua")]: false, + }), + node("p1", "Probe", "sst_probe", 0, { + fluid: "R134a", + tsat_c: 5.0, + [fixedFlagKey("tsat_c")]: true, + }), + ]; + const cfg = buildScenarioConfig(nodes, [ + { id: "pe", source: "e", target: "p1", sourceHandle: "outlet", targetHandle: "inlet" }, + ]); + expect(cfg.circuits[0].components[0].ua).toBe(8000); + expect(cfg.embeddings).toHaveLength(1); }); it("emits no embedding when Z_UA stays Fixed even with a Probe present", () => { diff --git a/apps/web/src/lib/configBuilder.ts b/apps/web/src/lib/configBuilder.ts index 7affe84..7fe99b0 100644 --- a/apps/web/src/lib/configBuilder.ts +++ b/apps/web/src/lib/configBuilder.ts @@ -310,8 +310,11 @@ export function buildScenarioConfig( const circuits = Array.from(circuitsMap.entries()) .sort(([a], [b]) => a - b) .map(([circuitId, cNodes]) => { + // Stable name order — React Flow node array order changes on drag/select. const components = cNodes .filter((n) => n.data.type !== CONTROL_NODE_TYPE) + .slice() + .sort((a, b) => a.data.name.localeCompare(b.data.name)) .map((n) => { const { type, name, params } = n.data; const canonicalParams = canonicalizeParams(type, params); @@ -394,9 +397,8 @@ export function buildScenarioConfig( * - `__fixed_*` Fixed checkbox flags * - measure-only targets (e.g. Probe `target`) — they become control setpoints * - emit Modelica-style `fix_pressure` / `fix_temperature` / `fix_mass_flow` - * - when Z_UA is Free (calibration), omit literal `ua` override — the CLI - * otherwise bakes `ua` into a fixed Calib factor and ignores live `z_ua`, - * which zeros ∂measure/∂z_ua and blows up Newton (singular J → bad Picard state) + * - BPHX only: when Z_UA is Free, omit absolute `ua` override (geometry = UA_nominal) + * - fill missing *required* params from catalogue defaults (e.g. FloodedEvaporator.ua) */ export function stripUiOnlyParams( type: string, @@ -417,12 +419,25 @@ export function stripUiOnlyParams( out[k] = v; } - // Free z_ua → live embedding owns UA scaling; drop absolute `ua` override. + // Free z_ua on BPHX: drop absolute `ua` override (geometry supplies UA_nominal). + // Legacy Condenser / Evaporator / Flooded*: keep `ua` — constructor nominal. const zUaMeta = meta?.params.find((p) => p.key === "z_ua" && p.actuatorFactor === "z_ua"); - if (zUaMeta && !isParamFixed(params, zUaMeta)) { + const bphxUaIsOverride = + type === "BphxCondenser" || type === "BphxEvaporator" || type === "BphxExchanger"; + if (bphxUaIsOverride && zUaMeta && !isParamFixed(params, zUaMeta)) { delete out.ua; } + // Required catalogue defaults (safety net when the canvas field was cleared). + for (const p of meta?.params ?? []) { + if (!p.required || p.default === undefined) continue; + if (bphxUaIsOverride && p.key === "ua") continue; + const cur = out[p.key]; + if (cur === undefined || cur === "") { + out[p.key] = p.default; + } + } + return applyExvFixSemantics(type, applyBoundaryFixSemantics(type, out, params), params); } @@ -596,13 +611,9 @@ export function buildModelEmbeddings( if (p.actuatorFactor && !fixed) { const n = typeof raw === "number" ? raw : Number(raw); - let initial = Number.isFinite(n) ? n : 1.0; - if ( - (p.actuatorFactor === "z_ua" || p.actuatorFactor === "z_dp") && - Math.abs(initial - 1.0) < 1e-12 - ) { - initial = 0.3; - } + // Keep the UI value as Newton start (do not rewrite 1.0 → 0.3 — that + // pushed Free Z into a fragile basin and amplified HashMap-order flakes). + const initial = Number.isFinite(n) ? n : 1.0; freeActs.push({ factor: p.actuatorFactor, initial, @@ -614,6 +625,7 @@ export function buildModelEmbeddings( } if (freeActs.length > 0) { + freeActs.sort((a, b) => a.factor.localeCompare(b.factor)); freeActsByComponent.set(node.data.name, { nodeName: node.data.name, node, acts: freeActs }); } } @@ -660,7 +672,10 @@ export function buildModelEmbeddings( }; const usedKeys = new Set(); - for (const { nodeName, acts } of freeActsByComponent.values()) { + const freeEntries = [...freeActsByComponent.values()].sort((a, b) => + a.nodeName.localeCompare(b.nodeName), + ); + for (const { nodeName, acts } of freeEntries) { for (const act of acts) { const probe = findProbeFor(nodeName, act.factor, usedKeys); if (probe) { @@ -683,6 +698,7 @@ export function buildModelEmbeddings( } } + embeddings.sort((a, b) => a.id.localeCompare(b.id)); return embeddings; } @@ -709,7 +725,7 @@ export function buildFixedFreeCalibrationControls( } /** Probe param key → semantic kind for pairing. */ -const PROBE_PARAM_KIND: Record = { +export const PROBE_PARAM_KIND: Record = { t_c: "T", tsat_c: "Tsat", p_bar: "P", @@ -779,8 +795,11 @@ function legacyProbeOutput(kind: string): string { } } -/** Probe physical kind → free factors it can calibrate. */ -const FACTOR_COMPATIBILITY: Record = { +/** + * Probe physical kind → free factors it can calibrate (Modelica Fixed ↔ Free). + * Shared by emit (`buildModelEmbeddings`) and UI aides (`calibAssist`). + */ +export const FACTOR_COMPATIBILITY: Record = { Tsat: ["z_ua"], Tsh: ["z_ua", "opening"], T: ["f_w"], diff --git a/crates/cli/src/config.rs b/crates/cli/src/config.rs index beabcee..12c8451 100644 --- a/crates/cli/src/config.rs +++ b/crates/cli/src/config.rs @@ -63,6 +63,12 @@ pub struct ScenarioConfig { /// (`output = value`). Distinct from [`Self::controls`] (no SaturatedController). #[serde(default)] pub embeddings: Vec, + /// Algebraic model equations: `lhs - rhs = 0` with `+ - * / min max exp ln`. + /// + /// Free unknown is either inferred from `lhs` when it is `component.z_*`, + /// or declared via [`EquationConfig::unknown`] (Probe target form). + #[serde(default)] + pub equations: Vec, /// Reusable subsystem templates (parameterized assemblies of components + /// internal edges, exposing a reduced external port set). Flattened into /// `circuits` at load time — the solver never sees the hierarchy. @@ -198,6 +204,56 @@ pub struct EmbeddingEquationConfig { pub value: f64, } +/// Algebraic model equation with a minimal expression language. +/// +/// ```text +/// lhs = rhs; // residual: lhs − rhs = 0 +/// ``` +/// +/// Supported ops: `+ - * /`, parentheses, unary `-`, `min`, `max`, `exp`, `ln`/`log`. +/// References: `component.field` or `"Component Name".field` +/// (`field` = `z_ua`, `T`, `Tsat`, `P`, …). +/// +/// **Form A** — free Z-factor on the left: +/// ```json +/// { "id": "eq1", "lhs": "evap.z_ua", "rhs": "min(1, exp(probe.T/300))", +/// "start": 0.3, "min": 0.05, "max": 2.0 } +/// ``` +/// +/// **Form B** — Probe target with explicit unknown (same role as embeddings[]): +/// ```json +/// { "id": "eq2", "lhs": "probe.Tsat", "rhs": "278.15", +/// "unknown": { "component": "evap", "factor": "z_ua", "start": 0.3, "min": 0.05, "max": 2.0 } } +/// ``` +#[derive(Debug, Clone, Serialize, Deserialize, JsonSchema)] +pub struct EquationConfig { + /// Unique id (also used as the residual / link id). + pub id: String, + /// Left-hand side expression. + pub lhs: String, + /// Right-hand side expression. + pub rhs: String, + /// Explicit free unknown (required when `lhs` is not itself a Z-factor ref). + #[serde(default)] + pub unknown: Option, + /// Start guess when the unknown is inferred from `lhs` (Form A). + #[serde(default = "default_actuator_initial")] + pub start: f64, + /// Lower bound when the unknown is inferred from `lhs` (Form A). + #[serde(default = "default_equation_min")] + pub min: f64, + /// Upper bound when the unknown is inferred from `lhs` (Form A). + #[serde(default = "default_equation_max")] + pub max: f64, +} + +fn default_equation_min() -> f64 { + 0.05 +} +fn default_equation_max() -> f64 { + 2.0 +} + /// A steady-state **system regulation** loop (EXV, injection, fan, …). /// /// Supports `SaturatedController`: saturated-PI with anti-windup, co-solved @@ -1399,6 +1455,7 @@ mod tests { "circuits", "controls", "embeddings", + "equations", "subsystems", "instances", "connections", @@ -1431,4 +1488,26 @@ mod tests { assert!((config.embeddings[0].equation.value - 278.15).abs() < 1e-9); assert!(config.controls.is_empty()); } + + #[test] + fn test_parse_equations() { + let json = r#"{ + "schema_version": "2", + "fluid": "R134a", + "equations": [{ + "id": "eq_evap_z_ua", + "lhs": "evap.z_ua", + "rhs": "min(1.0, exp(\"SST probe\".T / 300))", + "start": 0.3, + "min": 0.05, + "max": 2.0 + }] + }"#; + let config = ScenarioConfig::from_json(json).unwrap(); + assert_eq!(config.equations.len(), 1); + assert_eq!(config.equations[0].id, "eq_evap_z_ua"); + assert_eq!(config.equations[0].lhs, "evap.z_ua"); + assert!(config.equations[0].rhs.contains("min")); + assert!(config.equations[0].unknown.is_none()); + } } diff --git a/crates/cli/src/run.rs b/crates/cli/src/run.rs index 345b5d1..5b86e1e 100644 --- a/crates/cli/src/run.rs +++ b/crates/cli/src/run.rs @@ -646,23 +646,46 @@ fn execute_simulation( } } - // Free z_ua embedding: drop absolute `ua` overrides on the unknown's - // component. Otherwise `ua` freezes Calib and ∂Q/∂z_ua → 0 (singular J). - for emb in &config.embeddings { - if entropyk_core::normalize_factor_name(emb.unknown.factor.trim()) - == Some(entropyk_core::Z_UA) - { - if let Some(comp) = expanded_components - .iter_mut() - .find(|c| c.name == emb.unknown.component) + // Free z_ua on BPHX only: drop absolute `ua` override (geometry = UA_nominal). + // Legacy Condenser/Evaporator keep `ua` — it is the constructor nominal. + let mut free_z_ua: Vec<(String, String)> = config + .embeddings + .iter() + .filter(|e| { + entropyk_core::normalize_factor_name(e.unknown.factor.trim()) + == Some(entropyk_core::Z_UA) + }) + .map(|e| (e.unknown.component.clone(), e.id.clone())) + .collect(); + for eq in &config.equations { + if let Some(u) = eq.unknown.as_ref() { + if entropyk_core::normalize_factor_name(u.factor.trim()) == Some(entropyk_core::Z_UA) { - if comp.params.remove("ua").is_some() { - tracing::info!( - component = %comp.name, - embedding = %emb.id, - "Dropped 'ua' override — z_ua embedding owns UA scaling" - ); - } + free_z_ua.push((u.component.clone(), eq.id.clone())); + } + } else if let Ok(entropyk_solver::inverse::Expr::Ref { component, field }) = + entropyk_solver::inverse::parse_expr(&eq.lhs) + { + if entropyk_core::normalize_factor_name(&field) == Some(entropyk_core::Z_UA) { + free_z_ua.push((component, eq.id.clone())); + } + } + } + for (comp_name, src_id) in &free_z_ua { + if let Some(comp) = expanded_components + .iter_mut() + .find(|c| c.name == *comp_name) + { + let is_bphx = matches!( + comp.component_type.as_str(), + "BphxCondenser" | "BphxEvaporator" | "BphxExchanger" + ); + if is_bphx && comp.params.remove("ua").is_some() { + tracing::info!( + component = %comp.name, + source = %src_id, + "Dropped BPHX 'ua' override — free z_ua owns UA scaling" + ); } } } @@ -1151,6 +1174,34 @@ fn execute_simulation( } } + for eq in &config.equations { + match build_model_equation(eq) { + Ok((algebraic, bounded_var, unknown_id)) => { + if let Err(e) = system.add_algebraic_equation(algebraic) { + return control_error(format!( + "Failed to add algebraic equation '{}': {:?}", + eq.id, e + )); + } + if let Err(e) = system.add_bounded_variable(bounded_var) { + return control_error(format!( + "Failed to add unknown for equation '{}': {:?}", + eq.id, e + )); + } + if let Err(e) = system.link_constraint_to_control( + &entropyk_solver::inverse::ConstraintId::new(eq.id.clone()), + &unknown_id, + ) { + return control_error(format!("Failed to link equation '{}': {:?}", eq.id, e)); + } + } + Err(msg) => { + return control_error(format!("Invalid equation '{}': {}", eq.id, msg)); + } + } + } + for control in &config.controls { match build_saturated_control(control) { Ok((bounded_var, controller)) => { @@ -1527,6 +1578,20 @@ fn execute_simulation( } } } + for eq in &config.equations { + if let Ok((_, _, unknown_id)) = build_model_equation(eq) { + let start = eq + .unknown + .as_ref() + .map(|u| u.start) + .unwrap_or(eq.start); + if let Some(u_idx) = system.control_variable_state_index(&unknown_id) { + if u_idx < initial_state.len() { + initial_state[u_idx] = start; + } + } + } + } for control in &config.controls { let actuator_id = entropyk_solver::inverse::BoundedVariableId::new(saturated_actuator_id( &control.actuator.component, @@ -1597,6 +1662,7 @@ fn execute_simulation( .then(|| std::time::Duration::from_millis(config.solver.timeout_ms)); let needs_guarded_newton = !config.controls.is_empty() || !config.embeddings.is_empty() + || !config.equations.is_empty() || config.circuits.iter().any(|c| { c.enabled && c.components.iter().any(|comp| { @@ -3172,6 +3238,47 @@ fn bphx_calib_from_params( }) } +/// Calibration factors for legacy Condenser / Evaporator (`ua` is nominal, not override). +fn hx_calib_from_params( + params: &std::collections::HashMap, +) -> CliResult { + use entropyk_core::Calib; + let z_ua = params + .get("z_ua") + .or_else(|| params.get("Z_UA")) + .or_else(|| params.get("f_ua")) + .and_then(|v| v.as_f64()) + .unwrap_or(1.0); + if z_ua <= 0.0 { + return Err(CliError::Config(format!( + "z_ua must be > 0 (got {:.4})", + z_ua + ))); + } + let z_dp = params + .get("z_dp") + .or_else(|| params.get("Z_dpc")) + .or_else(|| params.get("f_dp")) + .and_then(|v| v.as_f64()) + .unwrap_or(1.0); + if z_dp <= 0.0 { + return Err(CliError::Config(format!( + "z_dp must be > 0 (got {:.4})", + z_dp + ))); + } + Ok(Calib { + z_flow: 1.0, + z_flow_eco: 1.0, + z_dp, + z_ua, + z_power: 1.0, + z_etav: 1.0, + f_w: 1.0, + calibration_source: None, + }) +} + /// Creates a pair of connected ports for components that need them (screw, MCHX, fan...). /// /// Ports are initialised at the given pressure and enthalpy. Both ports are connected @@ -3300,6 +3407,73 @@ fn build_model_embedding( Ok((constraint, bounded_var, unknown_id)) } +/// Algebraic model equation (`lhs − rhs = 0`) + free Z-factor unknown. +/// +/// Form A: `lhs` is `component.z_*` → unknown inferred, bounds from `start/min/max`. +/// Form B: explicit `unknown` block (Probe target: `lhs=probe.Tsat`, `rhs=315`). +fn build_model_equation( + eq: &crate::config::EquationConfig, +) -> Result< + ( + entropyk_solver::inverse::AlgebraicEquation, + entropyk_solver::inverse::BoundedVariable, + entropyk_solver::inverse::BoundedVariableId, + ), + String, +> { + use entropyk_solver::inverse::{ + parse_expr, AlgebraicEquation, BoundedVariable, BoundedVariableId, Expr, + }; + + let lhs = parse_expr(&eq.lhs).map_err(|e| format!("lhs: {e}"))?; + let rhs = parse_expr(&eq.rhs).map_err(|e| format!("rhs: {e}"))?; + + let (comp, factor, start, min, max) = if let Some(u) = eq.unknown.as_ref() { + let factor = u.factor.trim(); + if entropyk_core::normalize_factor_name(factor).is_none() { + return Err(format!( + "unknown Z-factor '{factor}' (expected z_ua, z_dp, z_flow, z_power, z_etav, f_w, …)" + )); + } + ( + u.component.clone(), + factor.to_string(), + u.start, + u.min, + u.max, + ) + } else if let Expr::Ref { component, field } = &lhs { + let canon = entropyk_core::normalize_factor_name(field).ok_or_else(|| { + format!( + "lhs '{component}.{field}' is not a Z-factor — set equations[].unknown \ + (e.g. Probe target form) or use lhs = component.z_ua" + ) + })?; + ( + component.clone(), + canon.to_string(), + eq.start, + eq.min, + eq.max, + ) + } else { + return Err( + "lhs must be component.z_* or provide equations[].unknown for the free Z-factor" + .into(), + ); + }; + + let unknown_id = BoundedVariableId::new(saturated_actuator_id(&comp, &factor)); + let bounded_var = BoundedVariable::with_component(unknown_id.clone(), &comp, start, min, max) + .map_err(|e| format!("invalid unknown bounds: {e:?}"))?; + let algebraic = AlgebraicEquation { + id: eq.id.clone(), + lhs, + rhs, + }; + Ok((algebraic, bounded_var, unknown_id)) +} + fn build_saturated_control( control: &crate::config::ControlConfig, ) -> Result< @@ -3596,7 +3770,17 @@ fn create_component( "FloodedEvaporator" => { use entropyk::FloodedEvaporator; - let ua = get_param_f64(params, "ua")?; + // UA is the nominal base (required for construction). Default matches + // the UI catalogue when the canvas field was left empty during calib. + let ua = params + .get("ua") + .and_then(|v| v.as_f64()) + .unwrap_or(8000.0); + if ua < 0.0 { + return Err(CliError::Config(format!( + "FloodedEvaporator: ua must be >= 0 (got {ua})" + ))); + } let target_quality = params .get("target_quality") .and_then(|v| v.as_f64()) @@ -3931,6 +4115,12 @@ fn create_component( cond = cond.with_flooded_head_pressure(target_k); } + // Fixed z_ua (or start hint): scale UA_eff = z_ua · UA. Live Free + // embeddings override via CalibIndices during solve. + if let Ok(calib) = hx_calib_from_params(params) { + cond.set_calib(calib); + } + Ok(Box::new(cond)) } @@ -4024,6 +4214,10 @@ fn create_component( evap = evap.with_regulated_superheat(); } + if let Ok(calib) = hx_calib_from_params(params) { + evap.set_calib(calib); + } + Ok(Box::new(evap)) } diff --git a/crates/components/src/heat_exchanger/condenser.rs b/crates/components/src/heat_exchanger/condenser.rs index 0fade68..a1dcb43 100644 --- a/crates/components/src/heat_exchanger/condenser.rs +++ b/crates/components/src/heat_exchanger/condenser.rs @@ -617,7 +617,19 @@ impl Condenser { FLOOD_LAMBDA_HI, FLOOD_LAMBDA_WIDTH, ); - jacobian.add_entry(row, act_idx, self.ua() * c.delta_t * c.e_exp * dlam); + let ua_nom = self.inner.ua_nominal() * self.live_z_ua(state); + jacobian.add_entry(row, act_idx, ua_nom * c.delta_t * c.e_exp * dlam); + } + // Live z_ua on secondary energy: r = ṁ·Δh − Q ⇒ ∂r/∂z = −∂Q/∂z + if let Some(z_ua_idx) = self.inner.calib_indices_ref().z_ua { + let flood_scale = if self.flood_ready() { + 1.0 - self.flooded_level(state) + } else { + 1.0 + }; + let d_q_dz = + self.inner.ua_nominal() * flood_scale * c.e_exp * c.delta_t; + jacobian.add_entry(row, z_ua_idx, -d_q_dz); } } None => { @@ -782,14 +794,27 @@ impl Condenser { } } - /// Effective conductance `UA_eff` [W/K] used by the coupled duty. With an - /// active flooding actuator this is `(1 − λ)·UA_nominal`; otherwise the - /// nominal `UA`. + /// Live calibration factor `z_ua` from Newton state when free, else calib param. + fn live_z_ua(&self, state: &StateSlice) -> f64 { + self.inner + .calib_indices_ref() + .z_ua + .and_then(|idx| state.get(idx).copied()) + .unwrap_or_else(|| self.calib().z_ua) + .max(0.0) + } + + /// Effective conductance `UA_eff` [W/K] used by the coupled duty. + /// + /// `UA_eff = UA_nominal · z_ua · (1 − λ)` when flooding is active, else + /// `UA_nominal · z_ua`. Reads live `state[z_ua]` when an embedding frees it + /// (otherwise ∂Q/∂z_ua = 0 → singular Jacobian). fn effective_ua(&self, state: &StateSlice) -> f64 { + let ua = self.inner.ua_nominal() * self.live_z_ua(state); if self.flood_ready() { - self.ua() * (1.0 - self.flooded_level(state)) + ua * (1.0 - self.flooded_level(state)) } else { - self.ua() + ua } } @@ -1592,7 +1617,7 @@ impl Component for Condenser { // ∂r1/∂λ = +UA_nom·(T_cond − T_sec,in)·e · dλ_eff/dλ. if self.flood_ready() { if let Some(lvl_idx) = self.fan_actuator_idx { - let ua_nom = self.ua(); + let ua_nom = self.inner.ua_nominal() * self.live_z_ua(state); let e = if c_sec > 1e-10 { (-ua_eff / c_sec).exp() } else { @@ -1608,6 +1633,24 @@ impl Component for Condenser { } } + // Live z_ua: UA = UA_nom·z_ua·(1−λ), ε = 1−e^(−UA/C), Q = ε·C·ΔT + // ⇒ ∂Q/∂z_ua = UA_nom·(1−λ)·e·ΔT ⇒ ∂r_energy/∂z_ua = −∂Q/∂z_ua + if let Some(z_ua_idx) = self.inner.calib_indices_ref().z_ua { + let e = if c_sec > 1e-10 { + (-ua_eff / c_sec).exp() + } else { + 0.0 + }; + let flood_scale = if self.flood_ready() { + 1.0 - self.flooded_level(state) + } else { + 1.0 + }; + let d_q_dz = + self.inner.ua_nominal() * flood_scale * e * (t_cond - t_sec_in); + jacobian.add_entry(row, z_ua_idx, -d_q_dz); + } + // r2 (emergent) = H_out − h_target(P_in): ∂/∂H_out = 1, // ∂/∂P_in = −dh_target/dP via central finite difference. if self.emergent_pressure { diff --git a/crates/components/src/heat_exchanger/evaporator.rs b/crates/components/src/heat_exchanger/evaporator.rs index 67df734..5bd7715 100644 --- a/crates/components/src/heat_exchanger/evaporator.rs +++ b/crates/components/src/heat_exchanger/evaporator.rs @@ -551,6 +551,18 @@ impl Evaporator { ); // ∂r/∂P_ref_in = ∂Q/∂P = −g·dT_evap/dP. jacobian.add_entry(row, c.ref_p_in_idx, -c.g * c.dtevap_dp); + // Live z_ua: secondary r = ṁ·Δh + Q ⇒ ∂r/∂z = +∂Q/∂z + if let Some(z_ua_idx) = self.inner.calib_indices_ref().z_ua { + let ua = self.live_ua(Some(state)); + let c_sec = state[m_in].abs() * c.cp_sec; + let e = if c_sec > 1e-10 && ua > 0.0 { + (-ua / c_sec).exp() + } else { + 0.0 + }; + let d_q_dz = self.inner.ua_nominal() * e * c.delta_t; + jacobian.add_entry(row, z_ua_idx, d_q_dz); + } } None => { jacobian.add_entry(row, h_out, 1.0); @@ -797,10 +809,22 @@ impl Evaporator { .map_err(|e| ComponentError::CalculationFailed(e.to_string())) } + /// Live UA = UA_nominal × z_ua (reads `state[calib_indices.z_ua]` when free). + fn live_ua(&self, state: Option<&StateSlice>) -> f64 { + let z = state + .and_then(|st| { + self.inner + .calib_indices_ref() + .z_ua + .and_then(|idx| st.get(idx).copied()) + }) + .unwrap_or_else(|| self.calib().z_ua); + self.inner.ua_nominal() * z.max(0.0) + } + /// Effectiveness for a phase-changing refrigerant (`C_min = C_sec`, `C_r → 0`): /// `ε = 1 − exp(−UA / C_sec)`. - fn effectiveness(&self, c_sec: f64) -> f64 { - let ua = self.ua(); + fn effectiveness(&self, c_sec: f64, ua: f64) -> f64 { if c_sec <= 1e-10 || ua <= 0.0 { return 0.0; } @@ -815,7 +839,8 @@ impl Evaporator { let c_sec = self.secondary_capacity_rate.unwrap_or(0.0); let t_sec_in = self.secondary_inlet_temp_k.unwrap_or(0.0); let t_evap = self.evap_temperature(p_in_pa)?; - let eps = self.effectiveness(c_sec); + let ua = self.live_ua(None); + let eps = self.effectiveness(c_sec, ua); Ok(eps * c_sec * (t_sec_in - t_evap)) } @@ -1053,7 +1078,8 @@ impl Component for Evaporator { // Live secondary stream: edge-driven in 4-port mode (Modelica). let (t_sec_in, c_sec) = self.live_secondary_stream(state)?; let t_evap = self.evap_temperature(p_in)?; - let eps = self.effectiveness(c_sec); + let ua = self.live_ua(Some(state)); + let eps = self.effectiveness(c_sec, ua); let q = eps * c_sec * (t_sec_in - t_evap); // r0: refrigerant pressure drop (tube MSH/Friedel + accel, or @@ -1266,7 +1292,8 @@ impl Component for Evaporator { // ∂r1/∂P_in = −∂Q/∂P_in = ε·C_sec·dT_evap/dP_in (T_sec,in constant), // dT_evap/dP via central finite difference. let (t_sec_in, c_sec) = self.live_secondary_stream(state)?; - let eps = self.effectiveness(c_sec); + let ua = self.live_ua(Some(state)); + let eps = self.effectiveness(c_sec, ua); let g = eps * c_sec; let t_evap = self.evap_temperature(p_in)?; let dp = p_in * 1e-4 + 100.0; @@ -1275,6 +1302,17 @@ impl Component for Evaporator { let dt_dp = (t_plus - t_minus) / (2.0 * dp); jacobian.add_entry(row, inlet_p_idx, g * dt_dp); + // Live z_ua: UA = UA_nom·z_ua, Q = ε·C·ΔT ⇒ ∂r_energy/∂z = −∂Q/∂z + if let Some(z_ua_idx) = self.inner.calib_indices_ref().z_ua { + let e = if c_sec > 1e-10 { + (-ua / c_sec).exp() + } else { + 0.0 + }; + let d_q_dz = self.inner.ua_nominal() * e * (t_sec_in - t_evap); + jacobian.add_entry(row, z_ua_idx, -d_q_dz); + } + // 4-port cross-derivatives of r1 to the secondary edge state: // ∂r1/∂h_sec,in = −∂Q/∂h_sec,in = −g·dT_sec/dh (exact 1/cp), // ∂r1/∂ṁ_sec = −∂Q/∂ṁ_sec = −g'(C_sec)·cp·(T_sec,in − T_evap). @@ -1282,7 +1320,6 @@ impl Component for Evaporator { let (m_s, p_s, h_s) = self.sec_in_idx.unwrap(); let cp_sec = self.sec_cp(state[p_s], state[h_s])?; let dt_dh = 1.0 / cp_sec; - let ua = self.ua(); let g_prime = if c_sec <= 1e-10 || ua <= 0.0 { 0.0 } else { @@ -2084,7 +2121,7 @@ mod tests { let c_sec = state[6] * cp_air; let t_air_in = (state[8] - 2_501_000.0 * w) / cp_air + 273.15; let t_evap = evap.evap_temperature(state[1]).unwrap(); - let eps = evap.effectiveness(c_sec); + let eps = evap.effectiveness(c_sec, evap.ua()); let q = eps * c_sec * (t_air_in - t_evap); assert!(q > 0.0, "evaporator must absorb heat: q={q}"); let expected = state[6] * (state[11] - state[8]) + q; diff --git a/crates/solver/src/inverse/bounded.rs b/crates/solver/src/inverse/bounded.rs index 485f0d7..9556e6c 100644 --- a/crates/solver/src/inverse/bounded.rs +++ b/crates/solver/src/inverse/bounded.rs @@ -42,7 +42,7 @@ use thiserror::Error; /// Type-safe identifier for a bounded control variable. /// /// Uses a string internally but provides type safety and clear intent. -#[derive(Debug, Clone, PartialEq, Eq, Hash)] +#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)] pub struct BoundedVariableId(String); impl BoundedVariableId { diff --git a/crates/solver/src/inverse/constraint.rs b/crates/solver/src/inverse/constraint.rs index 6a50f61..727eea0 100644 --- a/crates/solver/src/inverse/constraint.rs +++ b/crates/solver/src/inverse/constraint.rs @@ -16,7 +16,7 @@ use thiserror::Error; /// Type-safe identifier for a constraint. /// /// Uses a string internally but provides type safety and clear intent. -#[derive(Debug, Clone, PartialEq, Eq, Hash)] +#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)] pub struct ConstraintId(String); impl ConstraintId { diff --git a/crates/solver/src/inverse/embedding.rs b/crates/solver/src/inverse/embedding.rs index 08840c7..904a068 100644 --- a/crates/solver/src/inverse/embedding.rs +++ b/crates/solver/src/inverse/embedding.rs @@ -334,19 +334,24 @@ impl InverseControlConfig { self.control_to_constraint.get(bounded_variable_id) } - /// Returns an iterator over all constraint-to-control mappings. + /// Returns constraint→control mappings in **stable** order (sorted by constraint id). + /// + /// HashMap iteration is randomized per process; residual rows, Jacobian rows/columns, + /// and control state indices must share one deterministic order or Newton flakes. pub fn mappings(&self) -> impl Iterator { - self.constraint_to_control.iter() + let mut items: Vec<_> = self.constraint_to_control.iter().collect(); + items.sort_by(|a, b| a.0.cmp(b.0)); + items.into_iter() } - /// Returns an iterator over linked constraint IDs. + /// Linked constraint IDs in the same order as [`Self::mappings`]. pub fn linked_constraints(&self) -> impl Iterator { - self.constraint_to_control.keys() + self.mappings().map(|(c, _)| c) } - /// Returns an iterator over linked control variable IDs. + /// Linked control IDs in the same order as [`Self::mappings`] (column / state layout). pub fn linked_controls(&self) -> impl Iterator { - self.control_to_constraint.keys() + self.mappings().map(|(_, c)| c) } /// Checks if a constraint is linked. @@ -583,15 +588,18 @@ mod tests { fn test_inverse_control_config_mappings_iterator() { let mut config = InverseControlConfig::new(); + // Insert out of alphabetical order — iteration must still be sorted. config - .link(make_constraint_id("c1"), make_bounded_var_id("v1")) + .link(make_constraint_id("z_probe"), make_bounded_var_id("z_ua_cond")) .unwrap(); config - .link(make_constraint_id("c2"), make_bounded_var_id("v2")) + .link(make_constraint_id("a_probe"), make_bounded_var_id("z_ua_evap")) .unwrap(); - let mappings: Vec<_> = config.mappings().collect(); - assert_eq!(mappings.len(), 2); + let ids: Vec<&str> = config.mappings().map(|(c, _)| c.as_str()).collect(); + assert_eq!(ids, vec!["a_probe", "z_probe"]); + let controls: Vec<&str> = config.linked_controls().map(|c| c.as_str()).collect(); + assert_eq!(controls, vec!["z_ua_evap", "z_ua_cond"]); } #[test] diff --git a/crates/solver/src/inverse/expr.rs b/crates/solver/src/inverse/expr.rs new file mode 100644 index 0000000..105986e --- /dev/null +++ b/crates/solver/src/inverse/expr.rs @@ -0,0 +1,401 @@ +//! Minimal Modelica-style algebraic expression AST for model equations. +//! +//! Supports: `+ - * /`, parentheses, unary `-`, `min`, `max`, `exp`, `ln`/`log`, +//! numeric literals, and references `component.field` or `"Component Name".field`. + +use std::fmt; + +/// Parsed expression node. +#[derive(Debug, Clone, PartialEq)] +pub enum Expr { + Const(f64), + /// `component` + `field` (field is temperature, z_ua, …). + Ref { component: String, field: String }, + Neg(Box), + Add(Box, Box), + Sub(Box, Box), + Mul(Box, Box), + Div(Box, Box), + Min(Box, Box), + Max(Box, Box), + Exp(Box), + Ln(Box), +} + +impl Expr { + /// Evaluate with a resolver for `component.field` references. + pub fn eval(&self, resolve: &mut F) -> Result + where + F: FnMut(&str, &str) -> Result, + { + match self { + Expr::Const(v) => Ok(*v), + Expr::Ref { component, field } => resolve(component, field), + Expr::Neg(a) => Ok(-a.eval(resolve)?), + Expr::Add(a, b) => Ok(a.eval(resolve)? + b.eval(resolve)?), + Expr::Sub(a, b) => Ok(a.eval(resolve)? - b.eval(resolve)?), + Expr::Mul(a, b) => Ok(a.eval(resolve)? * b.eval(resolve)?), + Expr::Div(a, b) => { + let den = b.eval(resolve)?; + if den.abs() < 1e-30 { + return Err("division by zero in model equation".into()); + } + Ok(a.eval(resolve)? / den) + } + Expr::Min(a, b) => Ok(a.eval(resolve)?.min(b.eval(resolve)?)), + Expr::Max(a, b) => Ok(a.eval(resolve)?.max(b.eval(resolve)?)), + Expr::Exp(a) => { + let x = a.eval(resolve)?; + let y = x.exp(); + if y.is_finite() { + Ok(y) + } else { + Err(format!("exp({x}) overflow")) + } + } + Expr::Ln(a) => { + let x = a.eval(resolve)?; + if x <= 0.0 { + return Err(format!("ln({x}) domain error")); + } + Ok(x.ln()) + } + } + } + + /// Collect all `component.field` references. + pub fn collect_refs(&self, out: &mut Vec<(String, String)>) { + match self { + Expr::Const(_) => {} + Expr::Ref { component, field } => out.push((component.clone(), field.clone())), + Expr::Neg(a) | Expr::Exp(a) | Expr::Ln(a) => a.collect_refs(out), + Expr::Add(a, b) + | Expr::Sub(a, b) + | Expr::Mul(a, b) + | Expr::Div(a, b) + | Expr::Min(a, b) + | Expr::Max(a, b) => { + a.collect_refs(out); + b.collect_refs(out); + } + } + } +} + +/// Parse an expression string into an [`Expr`]. +pub fn parse_expr(input: &str) -> Result { + let mut p = Parser::new(input); + let e = p.parse_expr()?; + p.skip_ws(); + if p.pos < p.chars.len() { + return Err(format!( + "trailing junk at '{}'", + p.chars[p.pos..].iter().collect::() + )); + } + Ok(e) +} + +struct Parser { + chars: Vec, + pos: usize, +} + +impl Parser { + fn new(s: &str) -> Self { + Self { + chars: s.chars().collect(), + pos: 0, + } + } + + fn skip_ws(&mut self) { + while self.pos < self.chars.len() && self.chars[self.pos].is_whitespace() { + self.pos += 1; + } + } + + fn peek(&self) -> Option { + self.chars.get(self.pos).copied() + } + + fn bump(&mut self) -> Option { + let c = self.peek()?; + self.pos += 1; + Some(c) + } + + fn parse_expr(&mut self) -> Result { + let mut left = self.parse_term()?; + loop { + self.skip_ws(); + match self.peek() { + Some('+') => { + self.bump(); + let right = self.parse_term()?; + left = Expr::Add(Box::new(left), Box::new(right)); + } + Some('-') => { + self.bump(); + let right = self.parse_term()?; + left = Expr::Sub(Box::new(left), Box::new(right)); + } + _ => break, + } + } + Ok(left) + } + + fn parse_term(&mut self) -> Result { + let mut left = self.parse_unary()?; + loop { + self.skip_ws(); + match self.peek() { + Some('*') => { + self.bump(); + let right = self.parse_unary()?; + left = Expr::Mul(Box::new(left), Box::new(right)); + } + Some('/') => { + self.bump(); + let right = self.parse_unary()?; + left = Expr::Div(Box::new(left), Box::new(right)); + } + _ => break, + } + } + Ok(left) + } + + fn parse_unary(&mut self) -> Result { + self.skip_ws(); + if self.peek() == Some('-') { + self.bump(); + let e = self.parse_unary()?; + return Ok(Expr::Neg(Box::new(e))); + } + if self.peek() == Some('+') { + self.bump(); + return self.parse_unary(); + } + self.parse_primary() + } + + fn parse_primary(&mut self) -> Result { + self.skip_ws(); + match self.peek() { + Some('(') => { + self.bump(); + let e = self.parse_expr()?; + self.skip_ws(); + if self.bump() != Some(')') { + return Err("expected ')'".into()); + } + Ok(e) + } + Some('"') => self.parse_quoted_ref(), + Some(c) if c.is_ascii_digit() || c == '.' => self.parse_number(), + Some(c) if is_ident_start(c) => self.parse_ident_or_call(), + Some(c) => Err(format!("unexpected '{c}'")), + None => Err("unexpected end of expression".into()), + } + } + + fn parse_number(&mut self) -> Result { + let start = self.pos; + while let Some(c) = self.peek() { + if c.is_ascii_digit() || c == '.' || c == 'e' || c == 'E' { + self.bump(); + if (c == 'e' || c == 'E') && matches!(self.peek(), Some('+') | Some('-')) { + self.bump(); + } + } else { + break; + } + } + let s: String = self.chars[start..self.pos].iter().collect(); + s.parse::() + .map(Expr::Const) + .map_err(|_| format!("invalid number '{s}'")) + } + + fn parse_quoted_ref(&mut self) -> Result { + self.bump(); // " + let start = self.pos; + while let Some(c) = self.peek() { + if c == '"' { + break; + } + self.bump(); + } + let name: String = self.chars[start..self.pos].iter().collect(); + if self.bump() != Some('"') { + return Err("unterminated quoted component name".into()); + } + self.skip_ws(); + if self.bump() != Some('.') { + return Err("expected '.' after quoted component name".into()); + } + let field = self.parse_ident_string()?; + Ok(Expr::Ref { + component: name, + field, + }) + } + + fn parse_ident_string(&mut self) -> Result { + self.skip_ws(); + let start = self.pos; + if !self.peek().is_some_and(is_ident_start) { + return Err("expected identifier".into()); + } + self.bump(); + while self.peek().is_some_and(is_ident_cont) { + self.bump(); + } + Ok(self.chars[start..self.pos].iter().collect()) + } + + fn parse_ident_or_call(&mut self) -> Result { + let ident = self.parse_ident_string()?; + self.skip_ws(); + if self.peek() == Some('(') { + return self.parse_call(&ident); + } + // component.field[.more] — last segment is field, rest is component + let mut parts = vec![ident]; + while self.peek() == Some('.') { + self.bump(); + parts.push(self.parse_ident_string()?); + self.skip_ws(); + } + if parts.len() < 2 { + return Err(format!( + "bare identifier '{}' — expected component.field or function call", + parts[0] + )); + } + let field = parts.pop().unwrap(); + let component = parts.join("."); + Ok(Expr::Ref { component, field }) + } + + fn parse_call(&mut self, name: &str) -> Result { + self.bump(); // ( + self.skip_ws(); + let lower = name.to_ascii_lowercase(); + match lower.as_str() { + "min" | "max" => { + let a = self.parse_expr()?; + self.skip_ws(); + if self.bump() != Some(',') { + return Err(format!("expected ',' in {lower}(...)")); + } + let b = self.parse_expr()?; + self.skip_ws(); + if self.bump() != Some(')') { + return Err(format!("expected ')' after {lower}")); + } + Ok(if lower == "min" { + Expr::Min(Box::new(a), Box::new(b)) + } else { + Expr::Max(Box::new(a), Box::new(b)) + }) + } + "exp" => { + let a = self.parse_expr()?; + self.skip_ws(); + if self.bump() != Some(')') { + return Err("expected ')' after exp".into()); + } + Ok(Expr::Exp(Box::new(a))) + } + "ln" | "log" => { + let a = self.parse_expr()?; + self.skip_ws(); + if self.bump() != Some(')') { + return Err("expected ')' after ln".into()); + } + Ok(Expr::Ln(Box::new(a))) + } + other => Err(format!( + "unknown function '{other}' (supported: min, max, exp, ln)" + )), + } + } +} + +fn is_ident_start(c: char) -> bool { + c.is_ascii_alphabetic() || c == '_' +} + +fn is_ident_cont(c: char) -> bool { + c.is_ascii_alphanumeric() || c == '_' +} + +impl fmt::Display for Expr { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Expr::Const(v) => write!(f, "{v}"), + Expr::Ref { component, field } => write!(f, "{component}.{field}"), + Expr::Neg(a) => write!(f, "-({a})"), + Expr::Add(a, b) => write!(f, "({a} + {b})"), + Expr::Sub(a, b) => write!(f, "({a} - {b})"), + Expr::Mul(a, b) => write!(f, "({a} * {b})"), + Expr::Div(a, b) => write!(f, "({a} / {b})"), + Expr::Min(a, b) => write!(f, "min({a}, {b})"), + Expr::Max(a, b) => write!(f, "max({a}, {b})"), + Expr::Exp(a) => write!(f, "exp({a})"), + Expr::Ln(a) => write!(f, "ln({a})"), + } + } +} + +/// Algebraic model equation: `lhs - rhs = 0`. +#[derive(Debug, Clone)] +pub struct AlgebraicEquation { + pub id: String, + pub lhs: Expr, + pub rhs: Expr, +} + +impl AlgebraicEquation { + pub fn residual(&self, mut resolve: F) -> Result + where + F: FnMut(&str, &str) -> Result, + { + Ok(self.lhs.eval(&mut resolve)? - self.rhs.eval(&mut resolve)?) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn eval_const(e: &Expr) -> f64 { + e.eval(&mut |_, _| Err("no refs".into())).unwrap() + } + + #[test] + fn arithmetic_and_funcs() { + assert!((eval_const(&parse_expr("1 + 2 * 3").unwrap()) - 7.0).abs() < 1e-12); + assert!((eval_const(&parse_expr("min(2, 5)").unwrap()) - 2.0).abs() < 1e-12); + assert!((eval_const(&parse_expr("max(2, 5)").unwrap()) - 5.0).abs() < 1e-12); + assert!((eval_const(&parse_expr("exp(0)").unwrap()) - 1.0).abs() < 1e-12); + assert!((eval_const(&parse_expr("ln(exp(1))").unwrap()) - 1.0).abs() < 1e-9); + assert!((eval_const(&parse_expr("-(1+2)").unwrap()) + 3.0).abs() < 1e-12); + } + + #[test] + fn refs_and_quoted() { + let e = parse_expr(r#"evap.z_ua + "SST probe".T"#).unwrap(); + let v = e + .eval(&mut |c, f| match (c, f) { + ("evap", "z_ua") => Ok(0.1), + ("SST probe", "T") => Ok(280.0), + _ => Err("bad".into()), + }) + .unwrap(); + assert!((v - 280.1).abs() < 1e-12); + } +} diff --git a/crates/solver/src/inverse/mod.rs b/crates/solver/src/inverse/mod.rs index f0b7020..3d9f6b7 100644 --- a/crates/solver/src/inverse/mod.rs +++ b/crates/solver/src/inverse/mod.rs @@ -45,6 +45,7 @@ pub mod bounded; pub mod calibration; pub mod constraint; pub mod embedding; +pub mod expr; pub mod override_network; pub mod saturated_control; @@ -58,5 +59,6 @@ pub use calibration::{ }; pub use constraint::{ComponentOutput, Constraint, ConstraintError, ConstraintId}; pub use embedding::{ControlMapping, DoFError, InverseControlConfig}; +pub use expr::{parse_expr, AlgebraicEquation, Expr}; pub use override_network::{eval_error_signal, eval_error_weights, Combine, Objective}; pub use saturated_control::{SaturatedControlError, SaturatedController, Saturation}; diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 75a1bdb..3b07c74 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -16,7 +16,7 @@ use petgraph::algo; use petgraph::graph::{EdgeIndex, Graph, NodeIndex}; use petgraph::visit::EdgeRef; use petgraph::Directed; -use std::collections::HashMap; +use std::collections::{BTreeSet, HashMap}; use crate::coupling::{has_circular_dependencies, ThermalCoupling}; use crate::dof::{ @@ -208,6 +208,8 @@ pub struct System { thermal_couplings: Vec, /// Constraints for inverse control (output - target = 0) constraints: HashMap, + /// Modelica algebraic equations (`lhs - rhs = 0`), e.g. `evap.z_ua = f(...)`. + algebraic_equations: Vec, /// Bounded control variables for inverse control (with box constraints) bounded_variables: HashMap, /// Inverse control configuration (constraint → control variable mappings) @@ -246,6 +248,7 @@ impl System { node_to_circuit: HashMap::new(), thermal_couplings: Vec::new(), constraints: HashMap::new(), + algebraic_equations: Vec::new(), bounded_variables: HashMap::new(), inverse_control: InverseControlConfig::new(), saturated_controllers: Vec::new(), @@ -735,18 +738,20 @@ impl System { } } - if !self.constraints.is_empty() { + if !self.constraints.is_empty() || !self.algebraic_equations.is_empty() { match self.validate_inverse_control_dof() { Ok(()) => { tracing::debug!( - constraint_count = self.constraints.len(), + constraint_count = + self.constraints.len() + self.algebraic_equations.len(), control_count = self.inverse_control.mapping_count(), "Inverse control DoF validation passed" ); } Err(DoFError::UnderConstrainedSystem { .. }) => { tracing::warn!( - constraint_count = self.constraints.len(), + constraint_count = + self.constraints.len() + self.algebraic_equations.len(), control_count = self.inverse_control.mapping_count(), "Under-constrained inverse control system - solver may still converge" ); @@ -1240,6 +1245,21 @@ impl System { Ok(()) } + /// Adds a Modelica algebraic equation residual (`lhs − rhs = 0`). + pub fn add_algebraic_equation( + &mut self, + equation: crate::inverse::AlgebraicEquation, + ) -> Result<(), ConstraintError> { + let id = ConstraintId::new(equation.id.clone()); + if self.constraints.contains_key(&id) + || self.algebraic_equations.iter().any(|e| e.id == equation.id) + { + return Err(ConstraintError::DuplicateId { id }); + } + self.algebraic_equations.push(equation); + Ok(()) + } + /// Removes a constraint by ID. /// /// # Arguments @@ -1259,9 +1279,16 @@ impl System { self.constraints.len() } - /// Returns a reference to all constraints. + /// Number of Modelica algebraic equations (`lhs − rhs = 0`). + pub fn algebraic_equation_count(&self) -> usize { + self.algebraic_equations.len() + } + + /// Returns all constraints in **stable** order (sorted by constraint id). pub fn constraints(&self) -> impl Iterator { - self.constraints.values() + let mut items: Vec<_> = self.constraints.values().collect(); + items.sort_by(|a, b| a.id().cmp(b.id())); + items.into_iter() } /// Returns a reference to a specific constraint by ID. @@ -1317,8 +1344,12 @@ impl System { return Ok(0); } + // Must match Jacobian row order (`inverse_control.mappings()` / sorted ids). + let mut sorted: Vec<_> = self.constraints.values().collect(); + sorted.sort_by(|a, b| a.id().cmp(b.id())); + let mut count = 0; - for constraint in self.constraints.values() { + for constraint in sorted { let measured = match measured_values.get(constraint.id()).copied() { Some(v) => v, None => { @@ -2076,7 +2107,12 @@ impl System { constraint_id: &ConstraintId, bounded_variable_id: &BoundedVariableId, ) -> Result<(), DoFError> { - if !self.constraints.contains_key(constraint_id) { + let has_constraint = self.constraints.contains_key(constraint_id); + let has_algebraic = self + .algebraic_equations + .iter() + .any(|e| e.id == constraint_id.as_str()); + if !has_constraint && !has_algebraic { return Err(DoFError::ConstraintNotFound { constraint_id: constraint_id.clone(), }); @@ -2191,7 +2227,7 @@ impl System { pub fn validate_inverse_control_dof(&self) -> Result<(), DoFError> { let n_edge_unknowns = self.total_state_len; let n_controls = self.inverse_control.mapping_count(); - let n_constraints = self.constraints.len(); + let n_constraints = self.constraints.len() + self.algebraic_equations.len(); let n_unknowns = n_edge_unknowns + n_controls; let n_edge_eqs: usize = self @@ -2242,6 +2278,7 @@ impl System { .sum(); n_comp + self.constraints.len() + + self.algebraic_equations.len() + self.coupling_residual_count() + 2 * self.saturated_controllers.len() } @@ -2319,6 +2356,11 @@ impl System { name: id.to_string(), }); } + for eq in &self.algebraic_equations { + system_equations.push(EquationRole::ControlTracking { + name: eq.id.clone(), + }); + } for _ in 0..self.coupling_residual_count() { system_equations.push(EquationRole::CouplingDuty); } @@ -2701,6 +2743,7 @@ impl System { .map(|(_, c, _)| c.n_equations()) .sum(); total_eqs += self.constraints.len() + + self.algebraic_equations.len() + self.coupling_residual_count() + 2 * self.saturated_controllers.len(); @@ -2734,6 +2777,15 @@ impl System { .map_err(|e| ComponentError::CalculationFailed(e.to_string()))?; eq_offset += n_constraints; + // Modelica algebraic equations: lhs − rhs = 0 + for (i, eq) in self.algebraic_equations.iter().enumerate() { + let r = eq + .residual(|comp, field| self.resolve_expr_ref(comp, field, state)) + .map_err(ComponentError::CalculationFailed)?; + residuals[eq_offset + i] = r; + } + eq_offset += self.algebraic_equations.len(); + // Add couplings. Each coupling owns one unknown Q = state[coupling_state_index(i)]. // // * Physical mode (`hot_component`/`cold_component` set): the residual @@ -2836,12 +2888,19 @@ impl System { jacobian.add_entry(r, c, v); } + let algebraic_row_offset = row_offset + self.constraints.len(); + let algebraic_jac = self.compute_algebraic_equation_jacobian(state, algebraic_row_offset); + for (r, c, v) in algebraic_jac { + jacobian.add_entry(r, c, v); + } + // Thermal coupling rows: r_i = Q_i − η·duty_hot(state). // ∂r/∂Q = 1 always (also fixes the legacy stub's singular row); in // physical mode the plant coupling ∂r/∂col = −η·∂duty/∂col is formed by // central finite differences over the hot component's incident states // (same pattern as the saturated-controller measurement Jacobian). - let coupling_row_offset = row_offset + self.constraints.len(); + let coupling_row_offset = + row_offset + self.constraints.len() + self.algebraic_equations.len(); if !self.thermal_couplings.is_empty() { let eps = self.inverse_control.finite_diff_epsilon(); let mut state_mut = state.to_vec(); @@ -2882,8 +2941,10 @@ impl System { } } - let saturated_row_offset = - row_offset + self.constraints.len() + self.coupling_residual_count(); + let saturated_row_offset = row_offset + + self.constraints.len() + + self.algebraic_equations.len() + + self.coupling_residual_count(); let saturated_jac = self.compute_saturated_control_jacobian(state, saturated_row_offset); for (r, c, v) in saturated_jac { jacobian.add_entry(r, c, v); @@ -2891,6 +2952,147 @@ impl System { Ok(()) } + /// Resolve `component.field` for algebraic model equations. + /// + /// Fields may be Z-factors (`z_ua`, …) or measured outputs (`T`, `Tsat`, …). + pub fn resolve_expr_ref( + &self, + component: &str, + field: &str, + state: &StateSlice, + ) -> Result { + if let Some(canon) = entropyk_core::normalize_factor_name(field) { + let indices = self.calib_indices_by_name.get(component).ok_or_else(|| { + format!("unknown component '{component}' for Z-factor '{field}'") + })?; + let idx = match canon { + entropyk_core::Z_FLOW => indices.z_flow, + entropyk_core::Z_FLOW_ECO => indices.z_flow_eco, + entropyk_core::Z_DP => indices.z_dp, + entropyk_core::Z_UA => indices.z_ua, + entropyk_core::Z_POWER => indices.z_power, + entropyk_core::Z_ETAV => indices.z_etav, + entropyk_core::F_W => indices.f_w, + _ => None, + } + .ok_or_else(|| { + format!("Z-factor '{field}' is not a free unknown on component '{component}'") + })?; + return state.get(idx).copied().ok_or_else(|| { + format!("state index {idx} out of range for {component}.{field}") + }); + } + + let kind = Self::measured_output_from_field(field).ok_or_else(|| { + format!( + "unknown field '{field}' on '{component}' \ + (expected z_*, T/temperature, Tsat/SST/SDT, P/pressure, …)" + ) + })?; + let &node_idx = self.component_names.get(component).ok_or_else(|| { + format!("unknown component '{component}'") + })?; + let component_ref = self.graph.node_weight(node_idx).ok_or_else(|| { + format!("component '{component}' missing from graph") + })?; + component_ref + .measure_output(kind, state) + .filter(|v| v.is_finite()) + .ok_or_else(|| { + format!("measure '{field}' unavailable on component '{component}'") + }) + } + + fn measured_output_from_field(field: &str) -> Option { + use entropyk_components::MeasuredOutput as MO; + let n = field.trim().to_ascii_lowercase().replace(['_', '-'], ""); + Some(match n.as_str() { + "t" | "temperature" => MO::Temperature, + "tsat" | "sst" | "sdt" | "saturationtemperature" => MO::SaturationTemperature, + "p" | "pressure" => MO::Pressure, + "m" | "mdot" | "massflow" | "massflowrate" => MO::MassFlowRate, + "q" | "capacity" => MO::Capacity, + "duty" | "heattransferrate" | "heatrate" => MO::HeatTransferRate, + "superheat" | "tsh" | "sh" => MO::Superheat, + "subcooling" | "tsc" | "sc" => MO::Subcooling, + _ => return None, + }) + } + + fn compute_algebraic_equation_jacobian( + &self, + state: &StateSlice, + row_offset: usize, + ) -> Vec<(usize, usize, f64)> { + let mut entries = Vec::new(); + if self.algebraic_equations.is_empty() { + return entries; + } + let eps = self.inverse_control.finite_diff_epsilon(); + if state.len() < self.full_state_vector_len() { + tracing::error!( + state_len = state.len(), + required = self.full_state_vector_len(), + "compute_algebraic_equation_jacobian: state slice too short" + ); + return entries; + } + let mut state_mut = state.to_vec(); + + for (i, eq) in self.algebraic_equations.iter().enumerate() { + let row = row_offset + i; + let mut refs = Vec::new(); + eq.lhs.collect_refs(&mut refs); + eq.rhs.collect_refs(&mut refs); + + let mut cols = BTreeSet::new(); + for (comp, field) in &refs { + if let Some(canon) = entropyk_core::normalize_factor_name(field) { + if let Some(indices) = self.calib_indices_by_name.get(comp) { + let idx = match canon { + entropyk_core::Z_FLOW => indices.z_flow, + entropyk_core::Z_FLOW_ECO => indices.z_flow_eco, + entropyk_core::Z_DP => indices.z_dp, + entropyk_core::Z_UA => indices.z_ua, + entropyk_core::Z_POWER => indices.z_power, + entropyk_core::Z_ETAV => indices.z_etav, + entropyk_core::F_W => indices.f_w, + _ => None, + }; + if let Some(idx) = idx { + cols.insert(idx); + } + } + } + if let Some(&node_idx) = self.component_names.get(comp) { + for col in self.incident_state_indices_for_component(node_idx) { + cols.insert(col); + } + } + } + + for col in cols { + let orig = state_mut[col]; + state_mut[col] = orig + eps; + let r_plus = eq + .residual(|c, f| self.resolve_expr_ref(c, f, &state_mut)) + .unwrap_or(f64::NAN); + state_mut[col] = orig - eps; + let r_minus = eq + .residual(|c, f| self.resolve_expr_ref(c, f, &state_mut)) + .unwrap_or(f64::NAN); + state_mut[col] = orig; + if r_plus.is_finite() && r_minus.is_finite() { + let dr = (r_plus - r_minus) / (2.0 * eps); + if dr.abs() > 1e-10 { + entries.push((row, col, dr)); + } + } + } + } + entries + } + /// Tolerance for mass balance validation [kg/s]. /// /// This value (1e-9 kg/s) is tight enough to catch numerical issues