Buildings.Utilities.Math.Functions

Information

This package contains functions for commonly used mathematical operations. The functions are used in the blocks Buildings.Utilities.Math.

Package Content

Name	Description
bicubic	Bicubic function
biquadratic	Biquadratic function
cubicHermiteLinearExtrapolation	Interpolate using a cubic Hermite spline with linear extrapolation
polynomial	Polynomial function
powerLinearized	Power function that is linearized below a user-defined threshold
quadraticLinear	Function that is quadratic in first argument and linear in second argument
regNonZeroPower	Power function, regularized near zero, but nonzero value for x=0
smoothExponential	Once continuously differentiable approximation to exp(-|x|) in interval |x| < delta
smoothHeaviside	Once continuously differentiable approximation to the Heaviside function
smoothMax	Once continuously differentiable approximation to the maximum function
smoothMin	Once continuously differentiable approximation to the minimum function
smoothLimit	Once continuously differentiable approximation to the limit function
spliceFunction
splineDerivatives	Function to compute the derivatives for cubic hermite spline interpolation
trapezoidalIntegration	Integration using the trapezoidal rule
inverseXRegularized	Function that approximates 1/x by a twice continuously differentiable function
isMonotonic	Returns true if the argument is a monotonic sequence
Examples	Collection of models that illustrate model use and test models
BaseClasses	Package with base classes for Buildings.Utilities.Math.Functions

Buildings.Utilities.Math.Functions.bicubic

Information

y = a₁ + a₂ x₁ + a₃ x₁² + a₄ x₂ + a₅ x₂² + a₆ x₁ x₂ + a₇ x₁^3 + a₈ x₂^3 + a₉ x₁² x₂ + a₁0 x₁ x₂²

Inputs

Type	Name	Default	Description
Real	a[10]		Coefficients
Real	x1		Independent variable
Real	x2		Independent variable

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function bicubic "Bicubic function"
 input Real a[10] "Coefficients";
 input Real x1 "Independent variable";
 input Real x2 "Independent variable";
 output Real y "Result";
protected 
 Real x1Sq "= x1^2";
 Real x2Sq "= x2^2";
algorithm 
  x1Sq :=x1*x1;
  x2Sq :=x2*x2;
  y := a[1] + a[2] * x1 + a[3] * x1^2
            + a[4] * x2 + a[5] * x2^2
            + a[6] * x1 * x2
            + a[7] * x1Sq * x1
            + a[8] * x2Sq * x2
            + a[9] * x1Sq * x2
            + a[10] * x1 * x2Sq;

end bicubic;

Buildings.Utilities.Math.Functions.biquadratic

Information

y = a₁ + a₂ x₁ + a₃ x₁² + a₄ x₂ + a₅ x₂² + a₆ x₁ x₂

Inputs

Type	Name	Default	Description
Real	a[6]		Coefficients
Real	x1		Independent variable
Real	x2		Independent variable

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function biquadratic "Biquadratic function"
 input Real a[6] "Coefficients";
 input Real x1 "Independent variable";
 input Real x2 "Independent variable";
 output Real y "Result";
algorithm 
  y :=a[1] + x1*(a[2] + a[3]*x1) + x2*(a[4]+ a[5]*x2) + a[6]*x1*x2;

end biquadratic;

Buildings.Utilities.Math.Functions.cubicHermiteLinearExtrapolation

Information

For x₁ < x < x₂, this function interpolates using cubic hermite spline. For x outside this interval, the function linearly extrapolates.

For how to use this function, see Buildings.Utilities.Math.Functions.Examples.CubicHermite.

Inputs

Type	Name	Default	Description
Real	x		Abscissa value
Real	x1		Lower abscissa value
Real	x2		Upper abscissa value
Real	y1		Lower ordinate value
Real	y2		Upper ordinate value
Real	y1d		Lower gradient
Real	y2d		Upper gradient

Outputs

Type	Name	Description
Real	y	Interpolated ordinate value

Modelica definition

function cubicHermiteLinearExtrapolation 
  "Interpolate using a cubic Hermite spline with linear extrapolation"
  input Real x "Abscissa value";
  input Real x1 "Lower abscissa value";
  input Real x2 "Upper abscissa value";
  input Real y1 "Lower ordinate value";
  input Real y2 "Upper ordinate value";
  input Real y1d "Lower gradient";
  input Real y2d "Upper gradient";
  output Real y "Interpolated ordinate value";
algorithm 
  if (x > x1 and x < x2) then
    y:=Modelica.Fluid.Utilities.cubicHermite(
      x=x,
      x1=x1,
      x2=x2,
      y1=y1,
      y2=y2,
      y1d=y1d,
      y2d=y2d);
  elseif x <= x1 then
    // linear extrapolation
    y:=y1 + (x - x1)*y1d;
  else
    y:=y2 + (x - x2)*y2d;
  end if;
end cubicHermiteLinearExtrapolation;

Buildings.Utilities.Math.Functions.polynomial

Information

y = a₁ + a₂ x + a₃ x² + ...

Inputs

Type	Name	Default	Description
Real	x		Independent variable
Real	a[:]		Coefficients

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function polynomial "Polynomial function"
  annotation(derivative=Buildings.Utilities.Math.Functions.BaseClasses.der_polynomial);
 input Real x "Independent variable";
 input Real a[:] "Coefficients";
 output Real y "Result";
protected 
 parameter Integer n = size(a, 1)-1;
 Real xp[n+1] "Powers of x";
algorithm 
  xp[1] :=1;
  for i in 1:n loop
     xp[i+1] :=xp[i]*x;
  end for;
  y :=a*xp;
end polynomial;

Buildings.Utilities.Math.Functions.powerLinearized

Information

• the function is defined and monotone increasing for all x.
• dy/dx is bounded and continuous everywhere (for n < 1).

For x < x₀, this function replaces y=xⁿ by a linear function that is continuously differentiable everywhere.

A typical use of this function is to replace T = T4^(1/4) in a radiation balance to ensure that the function is defined everywhere. This can help solving the initialization problem when a solver may be far from a solution and hence T4 < 0.

See the package Examples for the graph.

Inputs

Type	Name	Default	Description
Real	x		Abscissa value
Real	n		Exponent
Real	x0		Abscissa value below which linearization occurs

Outputs

Type	Name	Description
Real	y	Function value

Modelica definition

function powerLinearized 
  "Power function that is linearized below a user-defined threshold"

 input Real x "Abscissa value";
 input Real n "Exponent";
 input Real x0 "Abscissa value below which linearization occurs";
 output Real y "Function value";
algorithm 
  if x > x0 then
   y := x^n;
  else
   y := x0^n * (1-n) + n * x0^(n-1) * x;
  end if;
end powerLinearized;

Buildings.Utilities.Math.Functions.quadraticLinear

Information

y = a₁ + a₂ x₁ + a₃ x₁² + (a₄ + a₅ x₁ + a₆ x₁²) x₂

Inputs

Type	Name	Default	Description
Real	a[6]		Coefficients
Real	x1		Independent variable for quadratic part
Real	x2		Independent variable for linear part

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function quadraticLinear 
  "Function that is quadratic in first argument and linear in second argument"
 input Real a[6] "Coefficients";
 input Real x1 "Independent variable for quadratic part";
 input Real x2 "Independent variable for linear part";
 output Real y "Result";
protected 
 Real x1Sq;
algorithm 
  x1Sq :=x1*x1;
  y :=a[1] + a[2]*x1 + a[3]*x1Sq + (a[4] + a[5]*x1 + a[6]*x1Sq)*x2;

end quadraticLinear;

Buildings.Utilities.Math.Functions.regNonZeroPower

Information

• y(0) is not equal to zero.
• dy/dx is bounded and continuous everywhere.

This function replaces y=|x|ⁿ in the interval -δ...+δ by a 4-th order polynomial that has the same function value and the first and second derivative at x=± δ.

A typical use of this function is to replace the function for the convective heat transfer coefficient for forced or free convection that is of the form h=c |dT|ⁿ for some constant c and exponent 0 ≤ n ≤ 1. By using this function, the original function that has an infinite derivative near zero and that takes on zero at the origin is replaced by a function with a bounded derivative and a non-zero value at the origin. Physically, the region -δ...+δ may be interpreted as the region where heat conduction dominates convection in the boundary layer.

Inputs

Type	Name	Default	Description
Real	x		Abscissa value
Real	n		Exponent
Real	delta	0.01	Abscissa value where transition occurs

Outputs

Type	Name	Description
Real	y	Function value

Modelica definition

function regNonZeroPower 
  "Power function, regularized near zero, but nonzero value for x=0"
  annotation(derivative=BaseClasses.der_regNonZeroPower);

 input Real x "Abscissa value";
 input Real n "Exponent";
 input Real delta = 0.01 "Abscissa value where transition occurs";
 output Real y "Function value";
protected 
  Real a1;
  Real a3;
  Real a5;
  Real delta2;
  Real x2;
  Real y_d "=y(delta)";
  Real yP_d "=dy(delta)/dx";
  Real yPP_d "=d^2y(delta)/dx^2";
algorithm 
  if abs(x) > delta then
   y := abs(x)^n;
  else
   delta2 :=delta*delta;
   x2 :=x*x;
   y_d :=delta^n;
   yP_d :=n*delta^(n - 1);
   yPP_d :=n*(n - 1)*delta^(n - 2);
   a1 := -(yP_d/delta - yPP_d)/delta2/8;
   a3 := (yPP_d - 12 * a1 * delta2)/2;
   a5 := (y_d - delta2 * (a3 + delta2 * a1));
   y := a5 + x2 * (a3 + x2 * a1);
   assert(a5>0, "Delta is too small for this exponent.");
  end if;
end regNonZeroPower;

Buildings.Utilities.Math.Functions.smoothExponential

Information

Function to provide a once continuously differentiable approximation to exp(- |x| ) in the interval |x| < δ for some positive δ

Inputs

Type	Name	Default	Description
Real	x		Input argument
Real	delta		Transition point where approximation occurs

Outputs

Type	Name	Description
Real	y	Output argument

Modelica definition

function smoothExponential 
  "Once continuously differentiable approximation to exp(-|x|) in interval |x| < delta"

  input Real x "Input argument";
  input Real delta "Transition point where approximation occurs";
  output Real y "Output argument";
protected 
  Real absX "Absolute value of x";
  Real a2 "Coefficient for approximating function";
  Real a3 "Coefficient for approximating function";
  Real e;
  Real d2;
  Real x2;
algorithm 
  absX :=abs(x);
  if absX > delta then
    y :=  exp(-absX);
  else
    d2 := delta*delta;
    e  := exp(-delta);
    a2 := (delta*e-4*(1-e))/2/d2;
    a3 := (e-1-a2*d2)/d2/d2;
    x2 := x*x;
    y  := 1+x2*(a2+x2*a3);
  end if;
end smoothExponential;

Buildings.Utilities.Math.Functions.smoothHeaviside

Information

Once Lipschitz continuously differentiable approximation to the Heaviside(.,.) function.

Inputs

Type	Name	Default	Description
Real	x		Argument
Real	delta		Parameter used for scaling

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function smoothHeaviside 
  "Once continuously differentiable approximation to the Heaviside function"
  input Real x "Argument";
  input Real delta "Parameter used for scaling";
  output Real y "Result";
algorithm 
 y := Buildings.Utilities.Math.Functions.spliceFunction(1, 0, x, delta);
end smoothHeaviside;

Buildings.Utilities.Math.Functions.smoothMax

Information

Once continuously differentiable approximation to the max(.,.) function.

Inputs

Type	Name	Default	Description
Real	x1		First argument
Real	x2		Second argument
Real	deltaX		Width of transition interval

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function smoothMax 
  "Once continuously differentiable approximation to the maximum function"
  input Real x1 "First argument";
  input Real x2 "Second argument";
  input Real deltaX "Width of transition interval";
  output Real y "Result";
algorithm 
  y := Buildings.Utilities.Math.Functions.spliceFunction(
         pos=x1, neg=x2, x=x1-x2, deltax=deltaX);
end smoothMax;

Buildings.Utilities.Math.Functions.smoothMin

Information

Once continuously differentiable approximation to the max(.,.) function.

Inputs

Type	Name	Default	Description
Real	x1		First argument
Real	x2		Second argument
Real	deltaX		Width of transition interval

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function smoothMin 
  "Once continuously differentiable approximation to the minimum function"
  input Real x1 "First argument";
  input Real x2 "Second argument";
  input Real deltaX "Width of transition interval";
  output Real y "Result";
algorithm 
  y := Buildings.Utilities.Math.Functions.spliceFunction(
       pos=x1, neg=x2, x=x2-x1, deltax=deltaX);
end smoothMin;

Buildings.Utilities.Math.Functions.smoothLimit

Information

Once continuously differentiable approximation to the limit(.,.) function. The output is bounded to be in [l, u].

Inputs

Type	Name	Default	Description
Real	x		Variable
Real	l		Low limit
Real	u		Upper limit
Real	deltaX		Width of transition interval

Outputs

Type	Name	Description
Real	y	Result

Modelica definition

function smoothLimit 
  "Once continuously differentiable approximation to the limit function"
  input Real x "Variable";
  input Real l "Low limit";
  input Real u "Upper limit";
  input Real deltaX "Width of transition interval";
  output Real y "Result";

protected 
  Real cor;
algorithm 
  cor :=deltaX/10;
  y := Buildings.Utilities.Math.Functions.smoothMax(x,l+deltaX,cor);
  y := Buildings.Utilities.Math.Functions.smoothMin(y,u-deltaX,cor);
end smoothLimit;

Buildings.Utilities.Math.Functions.spliceFunction

Information

Function to provide a once continuously differentialbe transition between to arguments.

The function is adapted from Modelica.Media.Air.MoistAir.Utilities.spliceFunction and provided here for easier accessability to model developers.

Inputs

Type	Name	Default	Description
Real	pos		Argument of x > 0
Real	neg		Argument of x < 0
Real	x		Independent value
Real	deltax		Half width of transition interval

Outputs

Type	Name	Description
Real	out	Smoothed value

Modelica definition

function spliceFunction
  annotation(derivative=BaseClasses.der_spliceFunction);
    input Real pos "Argument of x > 0";
    input Real neg "Argument of x < 0";
    input Real x "Independent value";
    input Real deltax "Half width of transition interval";
    output Real out "Smoothed value";
protected 
    Real scaledX1;
    Real y;
    constant Real asin1 = Modelica.Math.asin(1);
algorithm 
    scaledX1 := x/deltax;
    if scaledX1 <= -0.999999999 then
      out := neg;
    elseif scaledX1 >= 0.999999999 then
      out := pos;
    else
      y := (Modelica.Math.tanh(Modelica.Math.tan(scaledX1*asin1)) + 1)/2;
      out := pos*y + (1 - y)*neg;
    end if;

end spliceFunction;

Buildings.Utilities.Math.Functions.splineDerivatives

Information

This function computes the derivatives at the support points x_i that can be used as input for evaluating a cubic hermite spline. If ensureMonotonicity=true, then the support points y_i need to be monotone increasing (or increasing), and the computed derivatives d_i are such that the cubic hermite is monotone increasing (or decreasing). The algorithm to ensure monotonicity is based on the method described in Fritsch and Carlson (1980) for ρ = ρ₂.

This function is typically used with Buildings.Utilities.Math.Functions.cubicHermiteLinearExtrapolation which is used to evaluate the cubic spline. Because in many applications, the shape of the spline depends on parameters, this function has been implemented in such a way that all derivatives can be computed at once and then stored for use during the time stepping, in which the above function may be called.

References

F.N. Fritsch and R.E. Carlson, Monotone piecewise cubic interpolation. SIAM J. Numer. Anal., 17 (1980), pp. 238?246.

Inputs

Type	Name	Default	Description
Real	x[:]		Support point, strict monotone increasing
Real	y[size(x, 1)]		Function values at x
Boolean	ensureMonotonicity	isMonotonic(y, strict=false)	Set to true to ensure monotonicity of the cubic hermite

Outputs

Type	Name	Description
Real	d[size(x, 1)]	Derivative at the support points

Modelica definition

function splineDerivatives 
  "Function to compute the derivatives for cubic hermite spline interpolation"
  input Real x[:] "Support point, strict monotone increasing";
  input Real y[size(x, 1)] "Function values at x";
  input Boolean ensureMonotonicity=isMonotonic(y, strict=false) 
    "Set to true to ensure monotonicity of the cubic hermite";
  output Real d[size(x, 1)] "Derivative at the support points";
protected 
  Integer n=size(x, 1) "Number of data points";
  Real delta[n - 1] "Slope of secant line between data points";
  Real alpha "Coefficient to ensure monotonicity";
  Real beta "Coefficient to ensure monotonicity";
  Real tau "Coefficient to ensure monotonicity";

algorithm 
  if (n>1) then
    assert(x[1] < x[n], "x must be strictly increasing.
  Received x[1] = " + String(x[1]) + "
           x[" + String(n) + "] = " + String(x[n]));
  // Check data
    assert(isMonotonic(x, strict=true),
      "x-values must be strictly monontone increasing or decreasing.");
    if ensureMonotonicity then
      assert(isMonotonic(y, strict=false),
        "If ensureMonotonicity=true, y-values must be monontone increasing or decreasing.");
    end if;
  end if;

  // Compute derivatives at the support points
  if n == 1 then
    // only one data point
    d[1] :=0;
  elseif n == 2 then
    // linear function
    d[1] := (y[2] - y[1])/(x[2] - x[1]);
    d[2] := d[1];
  else
    // Slopes of the secant lines between i and i+1
    for i in 1:n - 1 loop
      delta[i] := (y[i + 1] - y[i])/(x[i + 1] - x[i]);
    end for;
    // Initial values for tangents at the support points.
    // End points use one-sided derivatives
    d[1] := delta[1];
    d[n] := delta[n - 1];

    for i in 2:n - 1 loop
      d[i] := (delta[i - 1] + delta[i])/2;
    end for;

  end if;
  // Ensure monotonicity
  if n > 2 and ensureMonotonicity then
    for i in 1:n - 1 loop
      if (abs(delta[i]) < Modelica.Constants.small) then
        d[i] := 0;
        d[i + 1] := 0;
      else
        alpha := d[i]/delta[i];
        beta := d[i + 1]/delta[i];
        // Constrain derivative to ensure monotonicity in this interval
        if (alpha^2 + beta^2) > 9 then
          tau := 3/(alpha^2 + beta^2)^(1/2);
          d[i] := delta[i]*alpha*tau;
          d[i + 1] := delta[i]*beta*tau;
        end if;
      end if;
    end for;
  end if;
end splineDerivatives;

Buildings.Utilities.Math.Functions.trapezoidalIntegration

Information

This function computes a definite integral using the trapezoidal rule.

Inputs

Type	Name	Default	Description
Integer	N		Number of integrand points
Real	f[:]		Integrands
Real	deltaX		Width of interval for Trapezoidal integration

Outputs

Type	Name	Description
Real	result	Result

Modelica definition

function trapezoidalIntegration 
  "Integration using the trapezoidal rule"
  input Integer N "Number of integrand points";
  input Real[:] f "Integrands";
  input Real deltaX "Width of interval for Trapezoidal integration";
  output Real result "Result";
algorithm 
  assert(N >= 2, "N must be no less than 2.");
  result := 0;
  for i in 1:N loop
    result := result + f[i];
  end for;

  result := 2*result;
  result := result - f[1] - f[N];
  result := result*deltaX/2;
end trapezoidalIntegration;

Buildings.Utilities.Math.Functions.inverseXRegularized

Information

Function that approximates y=1 ⁄ x inside the interval -δ ≤ x ≤ δ. The approximation is twice continuously differentiable with a bounded derivative on the whole real line.

See the package Examples for the graph.

Inputs

Type	Name	Default	Description
Real	x		Abscissa value
Real	delta		Abscissa value below which approximation occurs

Outputs

Type	Name	Description
Real	y	Function value

Modelica definition

function inverseXRegularized 
  "Function that approximates 1/x by a twice continuously differentiable function"
 input Real x "Abscissa value";
 input Real delta(min=0) "Abscissa value below which approximation occurs";
 output Real y "Function value";
protected 
 Real delta2 "Delta^2";
 Real x2_d2 "=x^2/delta^2";
algorithm 
  if (abs(x) > delta) then
    y := 1/x;
  else
    delta2 :=delta*delta;
    x2_d2  := x*x/delta2;
    y      := x/delta2 + x*abs(x/delta2/delta*(2 - x2_d2*(3 - x2_d2)));
  end if;

end inverseXRegularized;

Buildings.Utilities.Math.Functions.isMonotonic

Information

This function returns true if its argument is monotonic increasing or decreasing, and false otherwise. If strict=true, then strict monotonicity is tested, otherwise weak monotonicity is tested.

Inputs

Type	Name	Default	Description
Real	x[:]		Sequence to be tested
Boolean	strict	false	Set to true to test for strict monotonicity

Outputs

Type	Name	Description
Boolean	monotonic	True if x is monotonic increasing or decreasing

Modelica definition

function isMonotonic 
  "Returns true if the argument is a monotonic sequence"
  input Real x[:] "Sequence to be tested";
  input Boolean strict=false "Set to true to test for strict monotonicity";
  output Boolean monotonic "True if x is monotonic increasing or decreasing";
protected 
  Integer n=size(x, 1) "Number of data points";

algorithm 
  if n == 1 then
    monotonic := true;
  else
    monotonic := true;
    if strict then
      if (x[1] >= x[n]) then
        for i in 1:n - 1 loop
          if (not x[i] > x[i + 1]) then
            monotonic := false;
          end if;
        end for;
      else
        for i in 1:n - 1 loop
          if (not x[i] < x[i + 1]) then
            monotonic := false;
          end if;
        end for;
      end if;
    else
      // not strict
      if (x[1] >= x[n]) then
        for i in 1:n - 1 loop
          if (not x[i] >= x[i + 1]) then
            monotonic := false;
          end if;
        end for;
      else
        for i in 1:n - 1 loop
          if (not x[i] <= x[i + 1]) then
            monotonic := false;
          end if;
        end for;
      end if;
    end if;
    // strict
  end if;

end isMonotonic;