Skip to content

Commit f259145

Browse files
nachawatinachawati
committed
feat(core): improve zero-crossing bisection to 40 iterations, add chattering guard and time-event scheduling
1 parent 92f261e commit f259145

1 file changed

Lines changed: 107 additions & 7 deletions

File tree

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

Lines changed: 107 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,7 @@ export function generateFmuCSources(dae: ModelicaDAE, fmuResult: FmuResult, opti
6161
const fmi2FunctionsC = generateFmi2FunctionsC(id, nVars, nStates, nStringVars, dae, fmuResult);
6262

6363
// ── fmi3Functions.c ──
64-
const fmi3FunctionsC = generateFmi3FunctionsC(id);
64+
const fmi3FunctionsC = generateFmi3FunctionsC(id, dae);
6565

6666
// ── CMakeLists.txt ──
6767
const cmakeLists = generateCMakeLists(id);
@@ -260,6 +260,33 @@ function conditionToZeroCrossingC(condition: ModelicaExpression): string {
260260
return `(${exprToC(condition)} ? 1.0 : -1.0)`;
261261
}
262262

263+
/**
264+
* Extract time-event threshold from a when-condition.
265+
* Detects patterns like `time >= T` or `time > T` where T is a constant or expression.
266+
* Returns the C expression for the threshold, or null if the condition is not a time-event.
267+
*/
268+
function extractTimeEventThresholdC(condition: ModelicaExpression): string | null {
269+
if (!(condition instanceof ModelicaBinaryExpression)) return null;
270+
const op = condition.operator;
271+
// Pattern: time >= T or time > T
272+
if (
273+
(op === ModelicaBinaryOperator.GREATER_THAN_OR_EQUAL || op === ModelicaBinaryOperator.GREATER_THAN) &&
274+
condition.operand1 instanceof ModelicaNameExpression &&
275+
condition.operand1.name === "time"
276+
) {
277+
return exprToC(condition.operand2);
278+
}
279+
// Pattern: T <= time or T < time
280+
if (
281+
(op === ModelicaBinaryOperator.LESS_THAN_OR_EQUAL || op === ModelicaBinaryOperator.LESS_THAN) &&
282+
condition.operand2 instanceof ModelicaNameExpression &&
283+
condition.operand2.name === "time"
284+
) {
285+
return exprToC(condition.operand1);
286+
}
287+
return null;
288+
}
289+
263290
/**
264291
* Extract the assignment target variable name from the LHS of a simple equation.
265292
*/
@@ -642,6 +669,9 @@ function generateFmi2FunctionsC(
642669
lines.push(" fmi2Status asyncResult;");
643670
lines.push(" double asyncCurrentT;");
644671
lines.push(" double asyncTEnd;");
672+
lines.push(" /* Time-event scheduling fields */");
673+
lines.push(" int nextEventTimeDefined;");
674+
lines.push(" double nextEventTime;");
645675
lines.push("#ifdef _WIN32");
646676
lines.push(" HANDLE stepThread;");
647677
lines.push("#else");
@@ -920,11 +950,17 @@ function generateFmi2FunctionsC(
920950
lines.push(" double z0[N_EVENT_INDICATORS + 1];");
921951
lines.push(" double z1[N_EVENT_INDICATORS + 1];");
922952
}
953+
lines.push(" int eventCount = 0;");
923954

924955
lines.push(" while (t < tEnd - 1e-15 && !inst->cancelRequested) {");
925956
lines.push(" double step_h = h;");
926957
lines.push(" if (t + step_h > tEnd) step_h = tEnd - t;");
927958
lines.push("");
959+
// Check for scheduled time events and clamp step size
960+
lines.push(" if (inst->nextEventTimeDefined && inst->nextEventTime > t && inst->nextEventTime < t + step_h) {");
961+
lines.push(" step_h = inst->nextEventTime - t;");
962+
lines.push(" }");
963+
lines.push("");
928964

929965
for (let i = 0; i < stateRefs2.length; i++) {
930966
lines.push(` states0[${i}] = inst->model.vars[${stateRefs2[i]}];`);
@@ -947,9 +983,9 @@ function generateFmi2FunctionsC(
947983
lines.push(" if ((z0[i] > 0 && z1[i] <= 0) || (z0[i] <= 0 && z1[i] > 0)) { crossing = 1; break; }");
948984
lines.push(" }");
949985
lines.push(" if (crossing) {");
950-
lines.push(" /* Bisection root finding */");
986+
lines.push(" /* Bisection root finding (40 iterations ≈ 1e-12 precision) */");
951987
lines.push(" double h_left = 0, h_right = step_h;");
952-
lines.push(" for (int iter = 0; iter < 10; iter++) {");
988+
lines.push(" for (int iter = 0; iter < 40; iter++) {");
953989
lines.push(" double h_mid = 0.5 * (h_left + h_right);");
954990
lines.push(` take_rk4_step(&inst->model, t, h_mid, states0, states1);`);
955991
for (let i = 0; i < stateRefs2.length; i++) {
@@ -970,8 +1006,19 @@ function generateFmi2FunctionsC(
9701006
}
9711007
lines.push(" inst->model.time = t + step_h;");
9721008
lines.push("");
1009+
lines.push(" /* Chattering guard: prevent Zeno-type infinite event loops */");
1010+
lines.push(" eventCount++;");
1011+
lines.push(" if (eventCount > 100) {");
1012+
lines.push(" inst->asyncResult = fmi2Error;");
1013+
lines.push(" inst->asyncDone = 1;");
1014+
lines.push(" return; /* Abort: too many events in one communication step */");
1015+
lines.push(" }");
1016+
lines.push("");
9731017
lines.push(" fmi2EventInfo eventInfo;");
9741018
lines.push(" fmi2NewDiscreteStates((fmi2Component)inst, &eventInfo);");
1019+
lines.push(" /* Track scheduled time events from discrete state update */");
1020+
lines.push(" inst->nextEventTimeDefined = eventInfo.nextEventTimeDefined;");
1021+
lines.push(" if (eventInfo.nextEventTimeDefined) inst->nextEventTime = eventInfo.nextEventTime;");
9751022
lines.push(" }");
9761023
}
9771024

@@ -1199,6 +1246,33 @@ function generateFmi2FunctionsC(
11991246
}
12001247
}
12011248

1249+
// Time-event scheduling: scan when-conditions for "time >= T" patterns
1250+
// and set nextEventTimeDefined/nextEventTime to the nearest future event time
1251+
const timeThresholds: string[] = [];
1252+
for (const weq of whenEqs) {
1253+
const t = extractTimeEventThresholdC(weq.condition);
1254+
if (t) timeThresholds.push(t);
1255+
for (const clause of weq.elseWhenClauses) {
1256+
const t2 = extractTimeEventThresholdC(clause.condition);
1257+
if (t2) timeThresholds.push(t2);
1258+
}
1259+
}
1260+
if (timeThresholds.length > 0) {
1261+
lines.push("");
1262+
lines.push(" /* Scan for upcoming time events */");
1263+
lines.push(" {");
1264+
lines.push(" double t_current = inst->model.time;");
1265+
lines.push(" double t_next = 1e300;");
1266+
for (const threshold of timeThresholds) {
1267+
lines.push(` { double t_ev = ${threshold}; if (t_ev > t_current && t_ev < t_next) t_next = t_ev; }`);
1268+
}
1269+
lines.push(" if (t_next < 1e300) {");
1270+
lines.push(" info->nextEventTimeDefined = fmi2True;");
1271+
lines.push(" info->nextEventTime = t_next;");
1272+
lines.push(" }");
1273+
lines.push(" }");
1274+
}
1275+
12021276
lines.push(" return fmi2OK;");
12031277
lines.push("}");
12041278
lines.push("fmi2Status fmi2EnterContinuousTimeMode(fmi2Component c) { (void)c; return fmi2OK; }");
@@ -1561,7 +1635,7 @@ function collectReferencedNames(expr: unknown, names: Set<string>): void {
15611635
}
15621636
}
15631637
}
1564-
export function generateFmi3FunctionsC(id: string): string {
1638+
export function generateFmi3FunctionsC(id: string, dae: ModelicaDAE): string {
15651639
const lines: string[] = [];
15661640
lines.push("/* Auto-generated by ModelScript — FMI 3.0 API */");
15671641
lines.push(`#include "${id}_model.h"`);
@@ -1828,10 +1902,36 @@ export function generateFmi3FunctionsC(id: string): string {
18281902
lines.push(
18291903
" for (int i=0; i<n; i++) inst->model.states[i] = y0[i] + (h/6.0) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);",
18301904
);
1831-
lines.push(" if (eventHandlingNeeded) *eventHandlingNeeded = 0;");
1905+
1906+
// FMI 3.0 zero-crossing detection in fmi3DoStep
1907+
if (dae.eventIndicators.length > 0) {
1908+
lines.push("");
1909+
lines.push(" /* Zero-crossing detection */");
1910+
lines.push(" double z0_3[N_EVENT_INDICATORS+1], z1_3[N_EVENT_INDICATORS+1];");
1911+
lines.push(` /* Evaluate indicators at start (y0, t) */`);
1912+
lines.push(` for (int i=0; i<n; i++) inst->model.states[i] = y0[i];`);
1913+
lines.push(` inst->model.time = t;`);
1914+
lines.push(` ${id}_getEventIndicators(&inst->model, z0_3);`);
1915+
lines.push(` /* Restore end state */`);
1916+
lines.push(
1917+
` for (int i=0; i<n; i++) inst->model.states[i] = y0[i] + (h/6.0) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);`,
1918+
);
1919+
lines.push(` inst->model.time = t + h;`);
1920+
lines.push(` ${id}_getEventIndicators(&inst->model, z1_3);`);
1921+
lines.push(" int crossing3 = 0;");
1922+
lines.push(" for (int i=0; i<N_EVENT_INDICATORS; i++) {");
1923+
lines.push(" if ((z0_3[i] > 0 && z1_3[i] <= 0) || (z0_3[i] <= 0 && z1_3[i] > 0)) { crossing3 = 1; break; }");
1924+
lines.push(" }");
1925+
lines.push(" if (eventHandlingNeeded) *eventHandlingNeeded = crossing3;");
1926+
lines.push(" if (crossing3 && earlyReturn) *earlyReturn = 1;");
1927+
lines.push(" if (crossing3 && lastSuccessfulTime) *lastSuccessfulTime = t + h;");
1928+
} else {
1929+
lines.push(" if (eventHandlingNeeded) *eventHandlingNeeded = 0;");
1930+
lines.push(" if (earlyReturn) *earlyReturn = 0;");
1931+
lines.push(" if (lastSuccessfulTime) *lastSuccessfulTime = t + h;");
1932+
}
1933+
18321934
lines.push(" if (terminateSimulation) *terminateSimulation = 0;");
1833-
lines.push(" if (earlyReturn) *earlyReturn = 0;");
1834-
lines.push(" if (lastSuccessfulTime) *lastSuccessfulTime = t + h;");
18351935
lines.push(" return fmi3OK;");
18361936
lines.push("}");
18371937
lines.push("");

0 commit comments

Comments
 (0)