Skip to content

Commit 87ffd7c

Browse files
nachawatinachawati
committed
feat(core): symbolic differentiation engine, fmi2GetDirectionalDerivative jacobian codegen
1 parent 5fea06d commit 87ffd7c

3 files changed

Lines changed: 566 additions & 3 deletions

File tree

‎packages/core/src/compiler/modelica/fmi.ts‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -347,15 +347,15 @@ function generateModelDescriptionXml(
347347
if (opts.fmuType.modelExchange) {
348348
lines.push("");
349349
lines.push(
350-
` <ModelExchange modelIdentifier="${escapeXml(opts.modelIdentifier)}" canGetAndSetFMUstate="true" canSerializeFMUstate="true" />`,
350+
` <ModelExchange modelIdentifier="${escapeXml(opts.modelIdentifier)}" canGetAndSetFMUstate="true" canSerializeFMUstate="true" providesDirectionalDerivative="true" />`,
351351
);
352352
}
353353

354354
// CoSimulation element
355355
if (opts.fmuType.coSimulation) {
356356
lines.push("");
357357
lines.push(
358-
` <CoSimulation modelIdentifier="${escapeXml(opts.modelIdentifier)}" canGetAndSetFMUstate="true" canSerializeFMUstate="true" />`,
358+
` <CoSimulation modelIdentifier="${escapeXml(opts.modelIdentifier)}" canGetAndSetFMUstate="true" canSerializeFMUstate="true" providesDirectionalDerivative="true" />`,
359359
);
360360
}
361361

‎packages/core/src/compiler/modelica/fmu-codegen.ts‎

Lines changed: 84 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ import {
2929
ModelicaWhenEquation,
3030
} from "./dae.js";
3131
import type { FmuOptions, FmuResult } from "./fmi.js";
32+
import { differentiateExpr, simplifyExpr } from "./symbolic-diff.js";
3233
import { ModelicaBinaryOperator, ModelicaUnaryOperator, ModelicaVariability } from "./syntax.js";
3334

3435
/** Generated C source files. */
@@ -963,9 +964,91 @@ function generateFmi2FunctionsC(
963964
lines.push(" *state = (fmi2FMUstate)copy;");
964965
lines.push(" return fmi2OK;");
965966
lines.push("}");
967+
// ── fmi2GetDirectionalDerivative ──
968+
// Compute Δż = J · Δz where J is the Jacobian ∂f/∂x
969+
// J[i][j] = ∂(derivative[i])/∂(state[j]) — precomputed symbolically
966970
lines.push(
967-
"fmi2Status fmi2GetDirectionalDerivative(fmi2Component c, const fmi2ValueReference unknown[], size_t nUnknown, const fmi2ValueReference known[], size_t nKnown, const fmi2Real dvKnown[], fmi2Real dvUnknown[]) { (void)c; (void)unknown; (void)nUnknown; (void)known; (void)nKnown; (void)dvKnown; (void)dvUnknown; return fmi2Error; }",
971+
"fmi2Status fmi2GetDirectionalDerivative(fmi2Component c, const fmi2ValueReference unknown[], size_t nUnknown, const fmi2ValueReference known[], size_t nKnown, const fmi2Real dvKnown[], fmi2Real dvUnknown[]) {",
968972
);
973+
lines.push(" FMUInstance* inst = (FMUInstance*)c;");
974+
975+
// Extract derivative equations: der(x) = f(x, y, t)
976+
const derEquations: { stateName: string; rhs: ModelicaExpression }[] = [];
977+
for (const eq of dae.equations) {
978+
if (!("expression1" in eq && "expression2" in eq)) continue;
979+
const simpleEq = eq as { expression1: ModelicaExpression; expression2: ModelicaExpression };
980+
const lhsDer = extractDerName(simpleEq.expression1);
981+
const rhsDer = extractDerName(simpleEq.expression2);
982+
if (lhsDer) {
983+
derEquations.push({ stateName: lhsDer, rhs: simpleEq.expression2 });
984+
} else if (rhsDer) {
985+
derEquations.push({ stateName: rhsDer, rhs: simpleEq.expression1 });
986+
}
987+
}
988+
989+
// Build VR→derivative index and VR→state variable name mappings
990+
const derVRs = result.modelStructure.derivatives; // VRs of derivative variables
991+
const jacStateVRs: number[] = []; // VRs of state variables
992+
const stateNames: string[] = []; // Names of state variables
993+
for (const derEq of derEquations) {
994+
const sv = result.scalarVariables.find((v) => v.name === derEq.stateName);
995+
if (sv) {
996+
jacStateVRs.push(sv.valueReference);
997+
stateNames.push(derEq.stateName);
998+
}
999+
}
1000+
1001+
if (derEquations.length > 0 && jacStateVRs.length > 0) {
1002+
// Emit local variable aliases
1003+
lines.push(` double time = inst->model.time;`);
1004+
for (const sv of result.scalarVariables) {
1005+
if (sv.causality === "independent") continue;
1006+
const cName = varToC(sv.name);
1007+
lines.push(` double ${cName} = inst->model.vars[${sv.valueReference}];`);
1008+
}
1009+
lines.push("");
1010+
lines.push(" /* Zero output */");
1011+
lines.push(" for (size_t i = 0; i < nUnknown; i++) dvUnknown[i] = 0.0;");
1012+
lines.push("");
1013+
lines.push(" /* Accumulate Jacobian-vector product: dvUnknown[i] += J[i][j] * dvKnown[j] */");
1014+
lines.push(" for (size_t j = 0; j < nKnown; j++) {");
1015+
lines.push(" for (size_t i = 0; i < nUnknown; i++) {");
1016+
1017+
// For each (derivative, state) pair, emit the symbolic Jacobian entry
1018+
// Use switch on unknown VR, then switch on known VR
1019+
lines.push(" switch (unknown[i]) {");
1020+
for (let di = 0; di < derEquations.length; di++) {
1021+
const derVR = derVRs[di];
1022+
const derEq = derEquations[di];
1023+
if (derVR === undefined || !derEq) continue;
1024+
lines.push(` case ${derVR}: /* der(${derEq.stateName}) */`);
1025+
lines.push(" switch (known[j]) {");
1026+
for (let si = 0; si < stateNames.length; si++) {
1027+
const jacVR = jacStateVRs[si];
1028+
const stateName = stateNames[si];
1029+
if (jacVR === undefined || !stateName) continue;
1030+
// Symbolically differentiate rhs w.r.t. state variable
1031+
const jacobianEntry = simplifyExpr(differentiateExpr(derEq.rhs, stateName));
1032+
const jacobianC = exprToC(jacobianEntry);
1033+
lines.push(
1034+
` case ${jacVR}: dvUnknown[i] += (${jacobianC}) * dvKnown[j]; break; /* d/d(${stateName}) */`,
1035+
);
1036+
}
1037+
lines.push(" default: break;");
1038+
lines.push(" }");
1039+
lines.push(" break;");
1040+
}
1041+
lines.push(" default: break;");
1042+
lines.push(" }");
1043+
lines.push(" }");
1044+
lines.push(" }");
1045+
} else {
1046+
lines.push(" (void)inst; (void)unknown; (void)nUnknown; (void)known; (void)nKnown; (void)dvKnown;");
1047+
lines.push(" for (size_t i = 0; i < nUnknown; i++) dvUnknown[i] = 0.0;");
1048+
}
1049+
1050+
lines.push(" return fmi2OK;");
1051+
lines.push("}");
9691052
lines.push("");
9701053

9711054
return lines.join("\n");

0 commit comments

Comments
 (0)