Skip to content

Commit e9e8117

Browse files
nachawatinachawati
committed
feat(core): implement unified solver options, integrate kinsol wasm fallback for algebraic loops, wire solver configuration across fmu codegen and scripting api
1 parent 119ea2b commit e9e8117

12 files changed

Lines changed: 554 additions & 53 deletions

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

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717

1818
import { StaticTapeBuilder } from "./ad-codegen.js";
1919
import type { ModelicaExpression } from "./dae.js";
20+
import type { SolverOptions } from "./solver-options.js";
2021

2122
// ── Public interface ──
2223

@@ -34,6 +35,8 @@ export interface CoinorCodegenOptions {
3435
printLevel?: number;
3536
/** Use exact Hessian for IPOPT (default: true). */
3637
useExactHessian?: boolean;
38+
/** Generic solver options (overrides specific options). */
39+
solverOptions?: SolverOptions;
3740
}
3841

3942
/** Generated COIN-OR C files. */
@@ -96,9 +99,13 @@ export interface LpProblemDef {
9699
* Generate a CLP/CBC optimization driver from a LP/MILP problem definition.
97100
*/
98101
export function generateLpMainC(problem: LpProblemDef, options: CoinorCodegenOptions): CoinorCodegenResult {
99-
const solver = options.solver === "cbc" ? "cbc" : "clp";
100-
const tolerance = options.tolerance ?? 1e-8;
101-
const maxIter = options.maxIterations ?? 3000;
102+
const overrides = options.solverOptions;
103+
const resolvedSolver =
104+
overrides?.lpSolver === "clp" || overrides?.lpSolver === "cbc" ? overrides.lpSolver : options.solver;
105+
// If user asked for ipopt but we are in LP codegen, fallback to clp
106+
const solver = resolvedSolver === "cbc" ? "cbc" : "clp";
107+
const tolerance = overrides?.atol ?? options.tolerance ?? 1e-8;
108+
const maxIter = overrides?.maxNonlinearIterations ?? options.maxIterations ?? 3000;
102109
const printLevel = options.printLevel ?? 0;
103110

104111
const lines: string[] = [];
@@ -204,8 +211,9 @@ export function generateLpMainC(problem: LpProblemDef, options: CoinorCodegenOpt
204211
* Uses StaticTapeBuilder to emit AD-based objective/gradient/Hessian C code.
205212
*/
206213
export function generateNlpMainC(problem: NlpProblemDef, options: CoinorCodegenOptions): CoinorCodegenResult {
207-
const tolerance = options.tolerance ?? 1e-8;
208-
const maxIter = options.maxIterations ?? 3000;
214+
const overrides = options.solverOptions;
215+
const tolerance = overrides?.atol ?? options.tolerance ?? 1e-8;
216+
const maxIter = overrides?.maxNonlinearIterations ?? options.maxIterations ?? 3000;
209217
const printLevel = options.printLevel ?? 0;
210218

211219
const nVars = problem.variables.length;

‎packages/core/src/compiler/modelica/coinor-wasm.ts‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -385,7 +385,7 @@ export async function loadCoinorWasm(wasmUrl?: string): Promise<CoinorWasmSolver
385385
return cachedSolver;
386386
}
387387

388-
const url = wasmUrl ?? new URL("../../wasm/coinor.wasm", import.meta.url).href;
388+
const url = wasmUrl ?? new URL(/* webpackIgnore: true */ "../../wasm/coinor.wasm", import.meta.url).href;
389389
const jsUrl = url.replace(/\.wasm$/, ".js");
390390
const factory = await import(/* webpackIgnore: true */ jsUrl);
391391
const module = await new Promise<CoinorEmscriptenModule>((resolve) => {

‎packages/core/src/compiler/modelica/evaluate-optimize.ts‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ import {
1919
ModelicaStringLiteral,
2020
} from "./dae.js";
2121
import { ModelicaClassInstance } from "./model.js";
22+
import type { SolverOptions } from "./solver-options.js";
2223
import {
2324
ModelicaComponentReferenceSyntaxNode,
2425
type ModelicaFunctionCallSyntaxNode,
@@ -145,6 +146,21 @@ export function evaluateOptimize(
145146
const tolerance = getNamedArgNum("tolerance") ?? 1e-6;
146147
const maxIterations = getNamedArgNum("maxIterations") ?? 200;
147148

149+
// Parse solver options
150+
const solverOptions: SolverOptions = {};
151+
const _int = getNamedArgStr("integrator") as SolverOptions["integrator"];
152+
if (_int !== undefined) solverOptions.integrator = _int;
153+
const _nonlin = getNamedArgStr("nonlinear") as SolverOptions["nonlinear"];
154+
if (_nonlin !== undefined) solverOptions.nonlinear = _nonlin;
155+
const _lin = getNamedArgStr("linear") as SolverOptions["linear"];
156+
if (_lin !== undefined) solverOptions.linear = _lin;
157+
const _jac = getNamedArgStr("jacobian") as SolverOptions["jacobian"];
158+
if (_jac !== undefined) solverOptions.jacobian = _jac;
159+
const _opt = getNamedArgStr("optimizer") as SolverOptions["optimizer"];
160+
if (_opt !== undefined) solverOptions.optimizer = _opt;
161+
const _lp = getNamedArgStr("lpSolver") as SolverOptions["lpSolver"];
162+
if (_lp !== undefined) solverOptions.lpSolver = _lp;
163+
148164
// ── Step 4: Run optimization ──
149165
let result: {
150166
success: boolean;
@@ -166,6 +182,7 @@ export function evaluateOptimize(
166182
numIntervals,
167183
tolerance,
168184
maxIterations,
185+
solverOptions,
169186
});
170187
result = optimizer.optimize();
171188
} catch (e) {

‎packages/core/src/compiler/modelica/evaluate-simulate.ts‎

Lines changed: 31 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@ import {
2323
ModelicaStringLiteral,
2424
} from "./dae.js";
2525
import { ModelicaClassInstance } from "./model.js";
26+
import type { SolverOptions } from "./solver-options.js";
2627
import {
2728
ModelicaComponentReferenceSyntaxNode,
2829
type ModelicaFunctionCallSyntaxNode,
@@ -110,6 +111,16 @@ export function evaluateSimulate(
110111
return null;
111112
};
112113

114+
const getNamedArgStr = (name: string): string | null => {
115+
for (const na of namedArgs) {
116+
if (na.identifier?.text === name && na.argument?.expression) {
117+
const val = evaluateExpression(na.argument.expression, scope);
118+
if (val instanceof ModelicaStringLiteral) return val.value;
119+
}
120+
}
121+
return null;
122+
};
123+
113124
const exp = dae.experiment;
114125
const startTime = getNamedArg("startTime") ?? getPositionalArg(1) ?? exp.startTime ?? 0;
115126
const stopTime = getNamedArg("stopTime") ?? getPositionalArg(2) ?? exp.stopTime ?? 10;
@@ -121,13 +132,32 @@ export function evaluateSimulate(
121132
? outputIntervalArg
122133
: (exp.interval ?? (stopTime - startTime) / numberOfIntervals);
123134

135+
// Parse solver options
136+
const solverOptions: SolverOptions = {};
137+
const _int = getNamedArgStr("integrator") as SolverOptions["integrator"];
138+
if (_int !== undefined) solverOptions.integrator = _int;
139+
const _nonlin = getNamedArgStr("nonlinear") as SolverOptions["nonlinear"];
140+
if (_nonlin !== undefined) solverOptions.nonlinear = _nonlin;
141+
const _lin = getNamedArgStr("linear") as SolverOptions["linear"];
142+
if (_lin !== undefined) solverOptions.linear = _lin;
143+
const _jac = getNamedArgStr("jacobian") as SolverOptions["jacobian"];
144+
if (_jac !== undefined) solverOptions.jacobian = _jac;
145+
const _opt = getNamedArgStr("optimizer") as SolverOptions["optimizer"];
146+
if (_opt !== undefined) solverOptions.optimizer = _opt;
147+
const _lp = getNamedArgStr("lpSolver") as SolverOptions["lpSolver"];
148+
if (_lp !== undefined) solverOptions.lpSolver = _lp;
149+
124150
// ── Step 4: Run the simulation ──
125151
let result: { t: number[]; y: number[][]; states: string[] };
126152
let messages = "";
127153
try {
128154
const simulator = new deps.Simulator(dae);
129155
simulator.prepare();
130-
result = simulator.simulate(startTime, stopTime, step);
156+
result = simulator.simulate(startTime, stopTime, step, {
157+
atol: tolerance,
158+
rtol: tolerance,
159+
solverOptions,
160+
});
131161
} catch (e) {
132162
messages = e instanceof Error ? e.message : String(e);
133163
return buildResultRecord(startTime, stopTime, numberOfIntervals, tolerance, messages);

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ import {
2828
ModelicaSubscriptedExpression,
2929
ModelicaUnaryExpression,
3030
} from "./dae.js";
31+
import type { SolverOptions } from "./solver-options.js";
3132
import { ModelicaVariability } from "./syntax.js";
3233

3334
// ── Public interface ──
@@ -96,6 +97,8 @@ export interface FmuOptions {
9697
stepSize?: number | undefined;
9798
/** FMU type flags (default: both ME and CS). */
9899
fmuType?: FmuTypeFlags | undefined;
100+
/** Solver configuration. */
101+
solverOptions?: SolverOptions;
99102
}
100103

101104
/** Result of FMU generation. */

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

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -60,7 +60,8 @@ export function generateFmuCSources(dae: ModelicaDAE, fmuResult: FmuResult, opti
6060
const modelH = generateModelH(id, nVars, nStates, nStringVars, dae, fmuResult);
6161

6262
// ── model.c ──
63-
const modelC = generateModelC(id, dae, fmuResult) + "\n\n" + generateAlgebraicLoopSolvers(id, dae, fmuResult);
63+
const modelC =
64+
generateModelC(id, dae, fmuResult) + "\n\n" + generateAlgebraicLoopSolvers(id, dae, fmuResult, options);
6465

6566
// ── fmi2Functions.c ──
6667
const fmi2FunctionsC = generateFmi2FunctionsC(id, nVars, nStates, nStringVars, dae, fmuResult);
@@ -446,9 +447,15 @@ function generateModelH(
446447
return lines.join("\n");
447448
}
448449

449-
function generateAlgebraicLoopSolvers(id: string, dae: ModelicaDAE, result: FmuResult): string {
450+
function generateAlgebraicLoopSolvers(id: string, dae: ModelicaDAE, result: FmuResult, options: FmuOptions): string {
451+
const method = options.solverOptions?.jacobian ?? "ad-forward";
452+
450453
const lines: string[] = [];
451-
lines.push("/* Algebraic Loop Solver with Exact Analytical Jacobian (AD) */");
454+
if (method === "finite-difference") {
455+
lines.push("/* Algebraic Loop Solver with Finite-Difference Jacobian */");
456+
} else {
457+
lines.push("/* Algebraic Loop Solver with Exact Analytical Jacobian (AD) */");
458+
}
452459
lines.push(`#define LOG_ERROR(inst, msg) \\`);
453460
lines.push(` do { \\`);
454461
lines.push(` if ((inst)->logger) (inst)->logger((inst)->fmuInstance, "error", msg); \\`);

‎packages/core/src/compiler/modelica/init-solver.ts‎

Lines changed: 47 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,8 @@ import {
2121
ModelicaFunctionCallExpression,
2222
ModelicaNameExpression,
2323
} from "./dae.js";
24+
import { type SolverOptions } from "./solver-options.js";
25+
import { getCachedSundialsWasm } from "./sundials-wasm.js";
2426
import { ModelicaVariability } from "./syntax.js";
2527

2628
/** Result of initial equation solving. */
@@ -246,6 +248,7 @@ export function solveInitialEquations(
246248
startValues: Map<string, number>,
247249
parameters: Map<string, number>,
248250
startTime: number,
251+
solverOptions?: SolverOptions,
249252
): InitSolverResult {
250253
const result: InitSolverResult = {
251254
values: new Map(startValues),
@@ -368,13 +371,55 @@ export function solveInitialEquations(
368371
}
369372

370373
// Newton-Raphson iteration
371-
const maxIter = 50;
372-
const tol = 1e-10;
374+
const maxIter = solverOptions?.maxNonlinearIterations ?? 50;
375+
const tol = solverOptions?.atol ?? 1e-10;
373376
const nSolve = Math.min(nResiduals, nUnknowns); // Square system for Newton
374377

375378
// Initialize z from current env
376379
const z = unknownList.map((name) => env.get(name) ?? 0);
377380

381+
const useKinsol = solverOptions?.nonlinear === "kinsol" || solverOptions?.nonlinear === "hybrid";
382+
383+
if (useKinsol) {
384+
const solver = getCachedSundialsWasm();
385+
if (!solver) {
386+
throw new Error(
387+
"KINSOL solver requested but SUNDIALS WASM module is not loaded. Use simulateAsync() or loadSundialsWasm() first.",
388+
);
389+
}
390+
391+
const F = (zArr: number[]): number[] => {
392+
for (let i = 0; i < nUnknowns; i++) {
393+
const name = unknownList[i];
394+
if (name) env.set(name, zArr[i] ?? 0);
395+
}
396+
const res = new Array(nSolve);
397+
for (let row = 0; row < nSolve; row++) {
398+
const td = tapeData[row];
399+
if (!td) continue;
400+
const tArr = evaluateTapeForward(td.ops, env);
401+
res[row] = tArr[td.outputIndex] ?? 0;
402+
}
403+
return res;
404+
};
405+
406+
const kResult = solver.kinsol(F, z, { atol: tol, rtol: tol });
407+
if (kResult.converged || solverOptions?.nonlinear === "kinsol") {
408+
result.converged = kResult.converged;
409+
if (kResult.converged) {
410+
for (let i = 0; i < nUnknowns; i++) {
411+
const name = unknownList[i];
412+
if (name) result.values.set(name, kResult.solution[i] ?? 0);
413+
}
414+
}
415+
if (!kResult.converged && solverOptions?.nonlinear === "hybrid") {
416+
// Fall through to Newton-Raphson if hybrid and KINSOL failed
417+
} else {
418+
return result;
419+
}
420+
}
421+
}
422+
378423
for (let iter = 0; iter < maxIter; iter++) {
379424
result.iterations = iter + 1;
380425

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

Lines changed: 29 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
import type { ModelicaDAE } from "./dae.js";
1919
import { ModelicaIntegerLiteral, ModelicaRealLiteral } from "./dae.js";
2020
import { luFactor, luSolve, ModelicaSimulator } from "./simulator.js";
21+
import type { SolverOptions } from "./solver-options.js";
2122
import { ModelicaVariability } from "./syntax.js";
2223
import { Tape, type TapeNode } from "./tape.js";
2324

@@ -40,6 +41,10 @@ export interface OptimizationProblem {
4041
tolerance?: number;
4142
/** Maximum SQP iterations (default 200) */
4243
maxIterations?: number;
44+
/** Override parameters for the simulation */
45+
parameterOverrides?: Map<string, number>;
46+
/** Solver options for optimization and simulation */
47+
solverOptions?: SolverOptions;
4348
}
4449

4550
export interface OptimizationResult {
@@ -565,15 +570,38 @@ export class ModelicaOptimizer {
565570
evalJacobian,
566571
);
567572

573+
// Run one final simulation with the optimal controls to get the fine-grained state trajectories
574+
const optControls = new Map<string, number[]>();
575+
for (let j = 0; j < nControls; j++) {
576+
const name = controls[j]!;
577+
const uOpt = new Array<number>(nPoints);
578+
for (let k = 0; k < nPoints; k++) {
579+
uOpt[k] = result.z[k * varsPerPoint + nStates + j]!;
580+
}
581+
optControls.set(name, uOpt);
582+
}
583+
568584
// Extract results
569585
const stateTrajectories = new Map<string, number[]>();
570586
const controlTrajectories = new Map<string, number[]>();
571587

588+
const finalSimOpts = {
589+
parameterOverrides: new Map<string, number>(this.problem.parameterOverrides ?? []),
590+
...(this.problem.solverOptions ? { solverOptions: this.problem.solverOptions } : {}),
591+
};
592+
for (let k = 0; k < N; k++) {
593+
for (const name of controls) {
594+
finalSimOpts.parameterOverrides.set(name, optControls.get(name)![k]!);
595+
}
596+
}
597+
598+
const simResult = this.simulator.simulate(startTime, stopTime, dt, finalSimOpts);
599+
572600
for (let i = 0; i < nStates; i++) {
573601
const name = stateNames[i]!;
574602
const vals: number[] = [];
575603
for (let k = 0; k < nPoints; k++) {
576-
vals.push(result.z[k * varsPerPoint + i]!);
604+
vals.push(simResult.y[k]![i]!);
577605
}
578606
stateTrajectories.set(name, vals);
579607
}

0 commit comments

Comments
 (0)