19. Modelling Utilities

Like any derivative-based optimization solver, FORCESPRO works best if all functions defining the optimization problem are sufficiently smooth (i.e. at least continuously differentiable once). Both the MATLAB and the Python client of FORCESPRO come along with a couple of utility functions to assist the user with setting up such smooth problem formulations. This chapter provides details on those utility functions for modelling.

19.1. Interpolations 1D (e.g. splines)

If a given function is based on measurement data (or any other set of discrete data points), one can interpolate between those data points to yield a continuous function. FORCESPRO can either create such a function directly from the data points or allows you to provide a polynomial parameterization that can be used inside your symbolic problem formulation.

19.1.1. Polynomial Parameterization

A polynomial parameterization can be obtained by providing a vector containing M+1 break points (defining M interpolation intervals or pieces) as well as an array defining M sets of N+1 polynomial coefficients, each set defining a local polynomial of order N for each of those pieces.

Calling the line

MatlabPython

f = ForcesInterpolation(breaks, coefs);

will yield a symbolic representation of a polynomial in standard form

\[f(x) = \sum\limits_{j=0}^N c_{ij} x^j \quad \forall\,b_{i-1} \leq x \leq b_i \quad \forall\,i \in\{1,\ldots,M\}\,,\]

where b denotes the break points breaks and c denotes the coefficients coefs. f(x) can be a scalar or a K-dimensional function, i.e. coefs may be given for a multi-valued interpolation. For more details on how to pass those input parameters, we refer to the respective help function as the format differs slightly between the MATLAB and the Python client to follow domain-specific conventions.

In case your coefficients are defined relative to beginning of each piece, you can call

MatlabPython

f = ForcesInterpolation(breaks, coefs, 'pp');

to yield a symbolic representation of a polynomial in “piecewise polynomial” form

\[f(x) = \sum\limits_{j=0}^N c_{ij} (x-b_i)^j \quad \forall\,b_{i-1} \leq x \leq b_i \quad \forall\,i \in\{1,\ldots,M\}\,.\]

Note

In addition to providing fixed numerical values for break points and coefficients, you may also pass symbolic quantities for some or all of those! This will allow you to change the parameterization of your interpolation on the fly, e.g. by means of real-time parameters that are passed to the FORCESPRO solver.

The symbolic interpolation f can now be used inside your problem formulation by evaluating it, either at a fixed value or at any symbolic quantity, e.g.

MatlabPython

% assuming a state vector x and a control input u
y = f(x(1)) + u(1);

Important

Symbolic interpolations are currently only supported when using CasADi as AD tool.

19.1.2. Automatic Fit from Data

In case you do not want to specify break points and coefficients yourself, you can fit data points directly by calling:

MatlabPython

f = ForcesInterpolationFit(X, Y, method);

Here, X and Y are vectors (say, of dimension L) of data points to yield an interpolation that satisfies

\[f(X_i) = Y_i \quad \forall\,i \in\{1,\ldots,L\}\,.\]

The third argument method/kind specifies the method to be used to obtain that fit using built-in functionality of either MATLAB (see Table Table 19.1) or Python (see Table Table 19.2).

The symbolic interpolation f can be used the same way as described in section Section 19.1.1.

Table 19.1 Interpolation Method for MATLAB (see MATLAB’s interp1 for more details)

method	Description
'linear'	Piecewise linear
'nearest'	Piecewise constant, value from nearest data point
'next'	Piecewise constant, value from next data point
'previous'	Piecewise constant, value from previous data point
'spline' (default)	Piecewise cubic spline
'pchip'	Shape-preserving piecewise cubic spline

Table 19.2 Interpolation Method for Python (see SciPy’s interpolation class for more details)

kind	Description
'cubic' (default)	Piecewise cubic spline
'pchip'	Shape-preserving piecewise cubic spline

19.1.3. Application Example

A full example on how to use interpolations inside your problem formulation can be found here (MATLAB) or here (Python). Therein, both road limits are defined as splines and enforced as inequality constraints.

19.2. Interpolations 2D (B-splines)

In many industrial applications, 2D lookup tables/maps are commonly used. These have to be approximated by continuous function representations in order to use gradient-based solvers. Such lookup tables can be approximated using two-dimensional B-spline representations. A B-spline is a piecewise polynomial function.

In the following we provide several functionalities for fitting B-splines as well as capabilities for supplying coefficients of already existing B-splines directly. Importantly, all methods require the use of CasADi as the automatic differentiation tool. Furthermore, the symbolic variables have to be set to be CasADi MX expressions, which can be done by codeoptions.nlp.ad_expression_class = 'MX'.

19.2.1. B-Spline Representation

If you have already generated your own 2D spline, which can be represented as a general B-spline, you can directly provide the coefficients, knots and degrees. The coefficient matrix is matrix-indexed, meaning the x-inputs correspond to the rows and the y-inputs correspond to the columns. The knots positions knots_x and knots_y entail the internal knots as well as the additional ones on the boundaries/endpoints.

MatlabPython

f = ForcesInterpolation2D(knots_x, knots_y, coefs, degree_x, degree_y, m=1, casadiOptions=struct());

This function is a wrapper for casadi.Function.bspline, where m and casadiOptions are optional CasADi-function specific options. The returned f is a function handle taking as input the two inputs to the lookup table and returning the lookup table value: z=f(x, y).

19.2.2. Parametric B-Spline Representation

It is also possible to supply the coefficients online as runtime parameters. However, the knots and degrees still need to remain fixed, i.e. only the coefficient matrix can be changed online. The numel_coefs is the number of elements of the matrix-indexed coefficients matrix. The knots entail the internal as well as the additional ones on the boundaries/endpoints.

MatlabPython

f = ForcesInterpolation2DParametric(numel_coefs, knots_x, knots_y, degree_x, degree_y, m=1, casadiOptions=struct());

This function is a wrapper for casadi.bspline, where m and casadiOptions are optional CasADi-function specific options. The returned f is a function handle taking as inputs the two inputs to the lookup table as well as the coefficient matrix flattened in column-major order and returning the lookup table value: z=f(x, y, coefs).

19.2.3. Fitting B-Spline on Gridded Data

For fitting on grid data (i.e. rectangular mesh), we can directly use CasADi’s spline interpolation functionality.

MatlabPython

f = ForcesInterpolation2DGridded(xgrid, ygrid, Z, scaling=[1, 1, 1], casadiOptions=struct());

The xgrid and ygrid are the one-dimensional grid points in strictly ascending order and the Z is a two-dimensional matrix-indexed meshgrid of z values corresponding to xgrid and ygrid data sequences. Hence, Z has the size [length_x_grid, length_y_grid].

It additionally supports scaling of the spline inputs/outputs, which affects the numerical stability of the spline fitting routine. The scaling is automatically taken care of, so that the inputs/outputs to the resulting spline are still the physical quantities.

This function is a wrapper for casadi.interpolant, where casadiOptions are optional CasADi-function specific options. The returned f is a function handle taking as input the two inputs to the lookup table and returning the lookup table value: z=f(x, y).

The documentation for CasADi lookup tables can be found here: here (CasADi).

19.2.4. Fitting B-Spline on Arbitrary Data

For fitting on arbitrarily aligned data, we provide functionality which makes use of SciPy’s B-spline fitting routines: SmoothBivariateSpline and LSQBivariateSpline respectively. Hence, if utilizing these functions in MATLAB, a valid Python3 installation compatible with our MATLAB version needs to be present in our system. For a compatibility list of Python3 with different MATLAB versions, see here (Compatibility).

Note

Additionally, the Python3 package SciPy needs to be present in our system as well as the FORCESPRO client must be in our system’s Python path. For installing SciPy, follow the instructions in here (SciPy) and for instructions on how to add the FORCESPRO client to the Python path see here (Python Path).

The first method is the SmoothBivariateSpline of SciPy:

MatlabPython

f = ForcesInterpolationFit_SmoothBivariateSpline(x, y, z, w=ones(1, length(x)), bbox=[min(x), max(x), min(y), max(y)], kx=3, ky=3, s=length(w), rankTol=1e-16, scaling=[1, 1, 1], casadiOptions=struct());

All the arguments up to scaling are SciPy specific arguments, which can be found here: here (SmoothBivariateSpline).

It additionally supports scaling of the spline inputs/outputs, which affects the numerical stability of the spline fitting routine. The scaling is automatically taken care of, so that the inputs/outputs to the resulting spline are still the physical quantities.

The fitted coefficients of SciPy are provided to a wrapper for casadi.Function.bspline, where casadiOptions are optional CasADi-function specific options. The returned f is a function handle taking as input the two inputs to the lookup table and returning the lookup table value: z=f(x, y).

The second method is the LSQBivariateSpline of SciPy:

MatlabPython

f = ForcesInterpolationFit_LSQBivariateSpline(x, y, z, tx, ty, w=ones(1, length(x)), bbox=[min(x, tx), max(x, tx), min(y, ty), max(y, ty)], kx=3, ky=3, rankTol=1e-16, scaling=[1, 1, 1], casadiOptions=struct());

All the arguments up to scaling are SciPy specific arguments, which can be found here: here (LSQBivariateSpline).

It additionally supports scaling of the spline inputs/outputs, which affects the numerical stability of the spline fitting routine. The scaling is automatically taken care of, so that the inputs/outputs to the resulting spline are still the physical quantities.

The fitted coefficients of SciPy are provided to forcespro.modelling.Interpolation2D, which is a wrapper for casadi.Function.bspline, where casadiOptions are optional CasADi-function specific options. The returned f is a function handle taking as input the two inputs to the lookup table and returning the lookup table value: z=f(x, y).

19.2.5. Application Example

This example showcases how to generate splines from a 2D lookup table using the different provided functionalities.

Importantly, only the SMOOTH method was tuned properly, while LSQ and GRIDDED are tuned more roughly. Moreover, in this example, the coefficients provided to COEFF and PARAMETRIC are taken from the SMOOTH fitting solution.

The utilized 2D lookup represents the electric motor efficiency and is generated from figures in [JiaJibGor20_2D_modelling]. The 2D lookup table is stored as a matrix, where the rows correspond to the traction force, the columns to the velocity and each matrix entry represents the corresponding electric motor efficiency:

MatlabPython

spline_method = 'SMOOTH'; % 'SMOOTH', 'LSQ', 'COEFF', 'GRIDDED', 'PARAMETRIC'

vMax = 50;      % maximum velocity in the lookup table [m/s]
FtMax = 5e3;    % maximum traction force in the lookup table [N]

vSampleOrig = 0:vMax/10:vMax;
FtSampleOrig = 0:FtMax/5:FtMax;

vSampleOrig = [vSampleOrig(1), mean(vSampleOrig(1:2)), vSampleOrig(2:end)];
FtSampleOrig = [FtSampleOrig(1), mean(FtSampleOrig(1:2)), FtSampleOrig(2:end)];

% The data is taken from "Jia, Y.; Jibrin, R.; Goerges, D.: Energy-Optimal
% Adaptive Cruise Control for Electric Vehicles Based on Linear and
% Nonlinear Model Predictive Control. In: IEEE Transactions on Vehicular
% Technology, vol. 69, no. 12, pp. 14173-14187, Dec. 2020."
%     v [m/s] =        0,  2.5,    5,   10,   15,   20,   25,   30,   35,   40,   45,   50
etaSampleOrig = [
                    0.50, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50; ... % F [N] = 0
                    0.50, 0.68, 0.80, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.83, 0.81; ... % F [N] = 500
                    0.50, 0.78, 0.81, 0.88, 0.88, 0.88, 0.88, 0.88, 0.85, 0.83, 0.81, 0.80; ... % F [N] = 1000
                    0.50, 0.75, 0.81, 0.85, 0.88, 0.85, 0.83, 0.81, 0.78, 0.68, 0.65, 0.63; ... % F [N] = 2000
                    0.50, 0.68, 0.80, 0.83, 0.83, 0.81, 0.78, 0.68, 0.65, 0.63, 0.60, 0.57; ... % F [N] = 3000
                    0.50, 0.67, 0.78, 0.81, 0.80, 0.78, 0.65, 0.63, 0.60, 0.57, 0.55, 0.55; ... % F [N] = 4000
                    0.50, 0.67, 0.75, 0.78, 0.75, 0.65, 0.63, 0.60, 0.57, 0.55, 0.52, 0.52; ... % F [N] = 5000
                ];

% convert to 1-D sequences of data points
[vTemp, FtTemp] = meshgrid(vSampleOrig, FtSampleOrig);
vSample = reshape(vTemp.', 1, []);
FtSample = reshape(FtTemp.', 1, []);
etaSample = reshape(etaSampleOrig.', [1, numel(etaSampleOrig)]);

% spline
if strcmp(spline_method, 'SMOOTH')
    s = 0.025;
    w = ones(1, length(vSample));
    bbox = [min(vSample), max(vSample), min(FtSample), max(FtSample)];
    kx = 3;
    ky = 3;
    rankTol = 1e-8;
    scaling = [1e1, 1e3, 1];
    etaMotorSpline = ForcesInterpolationFit_SmoothBivariateSpline(vSample, FtSample, etaSample, w, bbox, kx, ky, s, rankTol, scaling);
elseif strcmp(spline_method, 'LSQ')
    tx = 30;
    ty = 3000;
    w = ones(1, length(vSample));
    bbox = [min([vSample, tx]), max([vSample, tx]), min([FtSample, ty]), max([FtSample, ty])];
    kx = 3;
    ky = 3;
    rankTol = 1e-8;
    scaling = [1e1, 1e3, 1];
    etaMotorSpline = ForcesInterpolationFit_LSQBivariateSpline(vSample, FtSample, etaSample, tx, ty, w, bbox, kx, ky, rankTol, scaling);
elseif strcmp(spline_method, 'COEFF')
    kx = 3;
    ky = 3;
    tx = [0, 0, 0, 0, 6.3993742001266, 24.1805094335686, 50, 50, 50, 50];
    ty = [0, 0, 0, 0, 760.320551795358, 1503.23831410745, 5000, 5000, 5000, 5000];
    % matrix-indexed coefficient matrix
    coeffs = [0.498727471092637   0.511037494098402   0.524901875945660   0.548511907792580   0.475607487652389   0.513135686437586
              0.498465143841756   0.654046675851965   0.783006359203673   0.713858506960411   0.676705972846789   0.681154627267599
              0.510442836827170   0.854158083414314   1.007125597517271   0.847628255554849   1.018658592375758   0.878826758082995
              0.495093430686428   0.735511289747710   0.857756122056489   0.863078750613390   0.548595365620131   0.497927393614425
              0.510240152112442   0.835399001169469   0.952683895958243   0.536982511482952   0.563982968586042   0.577416237725016
              0.501474795696226   0.773879473939183   0.878143062979889   0.444534437467682   0.615960539904494   0.508404545245928];

    etaMotorSpline = ForcesInterpolation2D(tx, ty, coeffs, kx, ky);
elseif strcmp(spline_method, 'GRIDDED')
    xgrid = vSampleOrig;
    ygrid = FtSampleOrig;
    Z = etaSampleOrig.'; % switch indexing: cartesian -> matrix
    scaling = [1e1, 1e3, 1];
    kx = 3;
    ky = 3;
    casadiOptions = struct('degree', [kx, ky]);

    etaMotorSpline = ForcesInterpolation2DGridded(xgrid, ygrid, Z, scaling, casadiOptions);
elseif strcmp(spline_method, 'PARAMETRIC')
    kx = 3;
    ky = 3;
    tx = [0, 0, 0, 0, 6.3993742001266, 24.1805094335686, 50, 50, 50, 50];
    ty = [0, 0, 0, 0, 760.320551795358, 1503.23831410745, 5000, 5000, 5000, 5000];
    % matrix-indexed coefficient matrix
    coeffs = [0.498727471092637   0.511037494098402   0.524901875945660   0.548511907792580   0.475607487652389   0.513135686437586
              0.498465143841756   0.654046675851965   0.783006359203673   0.713858506960411   0.676705972846789   0.681154627267599
              0.510442836827170   0.854158083414314   1.007125597517271   0.847628255554849   1.018658592375758   0.878826758082995
              0.495093430686428   0.735511289747710   0.857756122056489   0.863078750613390   0.548595365620131   0.497927393614425
              0.510240152112442   0.835399001169469   0.952683895958243   0.536982511482952   0.563982968586042   0.577416237725016
              0.501474795696226   0.773879473939183   0.878143062979889   0.444534437467682   0.615960539904494   0.508404545245928];

    numel_coeffs = numel(coeffs);
    etaMotorSpline = ForcesInterpolation2DParametric(numel_coeffs, tx, ty, kx, ky);
else
    error('Chosen spline_method is not valid');
end

% plotting
vMesh = 0:1:vMax;
FtMesh = 0:100:FtMax;
[vMesh, FtMesh] = meshgrid(vMesh, FtMesh);

if strcmp(spline_method, 'PARAMETRIC')
    etaMesh = full(etaMotorSpline(vMesh(:).', FtMesh(:).', coeffs(:))); % X/Y must be row-vectors
else
    etaMesh = full(etaMotorSpline(vMesh(:), FtMesh(:))); % sparse -> dense
end
etaMesh = reshape(etaMesh, size(vMesh, 1), size(vMesh, 2));

figure;
surf(vMesh, FtMesh, etaMesh);
xlabel("v [m/s]")
ylabel("Ft [N]")
title("Motor Efficiency")

figure;
contourf(vMesh, FtMesh, etaMesh, unique(etaSampleOrig)', 'ShowText',true)
xlabel("v [m/s]")
ylabel("Ft [N]")
title("Motor Efficiency")

[JiaJibGor20_2D_modelling]

Jia, Y.; Jibrin, R.; Goerges, D.: Energy-Optimal Adaptive Cruise Control for Electric Vehicles Based on Linear and Nonlinear Model Predictive Control. In: IEEE Transactions on Vehicular Technology, vol. 69, no. 12, pp. 14173-14187, Dec. 2020.

19.3. Smooth Approximations

There are a number of useful basic functions that are not differentiable everywhere. For some of them FORCESPRO provides a built-in smooth approximation and we plan to add more in an upcoming release.

19.3.1. Smooth Minimum

The minimum value of two scalars is not differentiable at the points where both values are identical. You can use the following smooth approximation instead:

MatlabPython

c = ForcesMin(a, b);

This function accepts an optional third argument to trade-off smoothness and approximation quality. The default value is set to 1e-8; higher values make the function smoother but less accurate.

19.3.2. Smooth Maximum

The maximum value of two scalars is not differentiable at the points where both values are identical. You can use the following smooth approximation instead:

MatlabPython

c = ForcesMax(a, b);

This function accepts an optional third argument to trade-off smoothness and approximation quality. The default value is set to 1e-8; higher values make the function smoother but less accurate.