为自定义Modelica插值函数添加显式导数以消除数值雅可比的问题

It looks like you misunderstood how the derivative annotation works in Modelica—you pointed it to a variable instead of a dedicated derivative function, which is why your changes didn’t resolve the warning. Let’s break down the correct approach to add an analytical derivative for your interpolate function and eliminate that numerical Jacobian.

What Went Wrong with Your Original Approach #

The derivative annotation in Modelica doesn’t reference a variable (like your dydx); it requires a separate function that calculates the derivative of your original function’s output with respect to its inputs. This function has a strict parameter format that must match the original function’s inputs, plus additional parameters for the original function’s output and the derivatives of each input.

Step-by-Step Fix for Linear Interpolation (Order=1) #

Since you’re using linear interpolation (the order=1 case), here’s how to implement the analytical derivative correctly:

1. Update the Original Interpolate Function’s Annotation #

First, modify your interpolate function to point its derivative annotation to a new derivative function (we’ll define this next):

function interpolate "Interpolate in a vector"
  extends Modelica.Icons.Function;
  input Real xVector[:] "Abscissa values (must be sorted)";
  input Real yVector[:] "Ordinate values (same length as xVector)";
  input Real x "Point to interpolate at";
  input Integer order=1 "Interpolation order (1=linear)";
  output Real y "Interpolated value";
  // Point to our upcoming derivative function
  annotation(derivative = interpolate_derivative);
end interpolate;

2. Define the Derivative Function #

Create a protected or sibling function interpolate_derivative that follows Modelica’s derivative function rules. For linear interpolation, the derivative of y with respect to x is just the slope of the segment containing x, scaled by the derivative of x itself. Here’s the implementation:

function interpolate_derivative
  // Match all inputs from the original interpolate function
  input Real xVector[:];
  input Real yVector[:];
  input Real x;
  input Integer order=1;
  // Add the original function's output (y)
  input Real y;
  // Add derivatives of each original input (most will be 0 for constants)
  input Real der_xVector[:];
  input Real der_yVector[:];
  input Real der_x;
  input Integer der_order;
  // Output the derivative of y (der_y)
  output Real der_y;
protected
  // Find the segment containing x (same logic as your original interpolate function)
  Integer idx = Modelica.Math.Vectors.findSegment(xVector, x);
  Real x1 = xVector[idx];
  Real x2 = xVector[idx+1];
  Real y1 = yVector[idx];
  Real y2 = yVector[idx+1];
  Real slope;
algorithm
  if order == 1 then
    // Calculate slope (handle near-zero segment length to avoid division by zero)
    if abs(x2 - x1) > Modelica.Constants.eps then
      slope = (y2 - y1)/(x2 - x1);
    else
      slope = sign(y2 - y1)*Modelica.Constants.inf;
    end if;
    // For linear interpolation, der_y = slope * der_x
    // (If xVector/yVector were time-varying, you'd add terms for der_xVector/der_yVector here)
    der_y = slope * der_x;
  else
    // Add handling for higher orders if you need them; for now, default to 0
    der_y = 0;
  end if;
end interpolate_derivative;

3. Verify the Fix #

After implementing these changes, re-translate your model in Dymola. The Cannot find derivative of function warning should disappear, and Dymola will use your analytical derivative instead of generating a numerical Jacobian.

Key Notes #

• Parameter Matching: The derivative function’s input order must exactly match the original function’s inputs, followed by the original function’s output, then the derivatives of each original input.
• Constant Inputs: If xVector and yVector are constants (which they often are for interpolation tables), their derivatives (der_xVector, der_yVector) will be 0, so they don’t affect the result.
• Higher Orders: If you ever need to support higher-order interpolation (like cubic splines), you’ll need to extend the derivative function to calculate the appropriate analytical derivative for that order.

内容的提问来源于stack exchange，提问作者Justin Kauffman