Programmatic Surface Fitting :: Programmatic Curve and Surface Fitting (Curve Fitting Toolbox™)

Programmatic Surface Fitting

On this page…
Surface Fitting Objects and Methods Automotive Fuel Efficiency Programmatic Surface Fitting Example Biopharmaceutical Drug Interaction Programmatic Surface Fitting Example

On this page…

Automotive Fuel Efficiency Programmatic Surface Fitting Example

Biopharmaceutical Drug Interaction Programmatic Surface Fitting Example

Surface Fitting Objects and Methods

This section describes how to use Curve Fitting Toolbox functions from the command-line or to write programs for surface fitting applications.

Overview

One way to quickly assemble code for surface fits and plots into useful programs is to generate an M-file from a session in Surface Fitting Tool. In this way, you can transform your interactive analysis of a single data set into a reusable function for command-line analysis or for batch processing of multiple data sets. You can use the generated M-file without modification, or edit and customize the code as needed. See Generating M-Files from the Surface Fitting Tool.

Surface Fitting Objects and Methods

The surface fit object (sfit) stores the results from a surface fitting operation, making it easy to plot and analyze fits at the command line.

Like cfit objects, sfit objects are a subclass of fittype objects, so they inherit all the same methods of fittype. For a list of available methods, see Fit Type Methods.

sfit objects also provide methods exclusively for sfit objects, listed in the table here: Surface Fit Methods.

See Function Reference for information on all Curve Fitting Toolbox functions, classes, and methods.

Automotive Fuel Efficiency Programmatic Surface Fitting Example

Curve Fitting Toolbox software provides some example data generated from a GTPOWER predictive combustion engine model. The model emulates a naturally aspirated spark ignition, 2-liter, inline 4-cylinder engine. You can use surface fitting methods to fit a response surface to this data to investigate fuel efficiency.

The data set includes the following variables you need to model response surfaces:

Speed is in revolutions per minute (RPM) units.
Load is the normalized cylinder air mass (the ratio of cylinder aircharge to maximum naturally aspirated cylinder aircharge at standard temperature and pressure).
BSFC is the brake-specific fuel consumption in g/Kwh; that is, the energy flow in divided by mechanical power out (fuel efficiency).

The aim is to model a response surface to find the minimum BSFC as a function of speed and load. You can use this surface as a table, included as part of a hybrid vehicle optimization algorithm combining the use of a motor and your engine. To operate the engine as fuel efficiently as possible, the table must operate the engine near the bottom of the BSFC bowl.

Load and Preprocess Data

Follow these steps to load and process the data:

Load the data from the XLS spreadsheet. Use the 'basic' command option for non- Windows^® platforms.
Create a variable n that has all the numeric data in one array.
```
n = xlsread( 'Engine_Data_SI_NA_2L_I4.xls', 'SI NA 2L I4',...
 '', 'basic' );
```

Extract from the variable n the columns of interest:

SPEED = n(:,2);
LOAD_CMD = n(:,3);
LOAD = n(:,8);
BSFC = n(:,22);

Process the data before fitting, to pick out the min(BSFC) values from each sweep. The data points are organized insweeps on speed/load.

Get a list of the speed/load sites:

SL = unique( [SPEED, LOAD_CMD], 'rows' );
nRuns = size( SL, 1 );

For each speed/load site, find the data at the site and extract the actual measured load and the minimum BSFC.

minBSFC  = zeros( nRuns, 1 );
Load     = zeros( nRuns, 1 );
Speed    = zeros( nRuns, 1 );
for i = 1:nRuns
    idx = SPEED == SL(i,1) & LOAD_CMD == SL(i,2);

    minBSFC(i) = min(  BSFC(idx) );
    Load(i)    = mean( LOAD(idx) );
    Speed(i)   = mean( SPEED(idx) );
end

Fit and Plot Surfaces of Fuel Efficiency

Follow these steps to fit and plot some surfaces:

Fit a surface to the preprocessed data.

f1 = fit( [Speed, Load], minBSFC, 'Lowess', 'Normalize', 'on' )

This command results in the following output:

     Locally weighted smoothing linear regression:
       f1(x,y) = lowess (linear) smoothing regression
       computed from p
       where x is normalized by mean 3407 and std 1214
       and where y is normalized by mean 0.5173 and std 0.1766
     Coefficients:
       p = coefficient structure

Plot your fit:

plot( f1, [Speed, Load], minBSFC );
xlabel( 'Speed [RPM]' );
ylabel( 'Load [%]' );
zlabel( 'Minimum BSFC [g/Kwh]' );

Review the resulting plot:

There are points where BSFC is negative because this data is generated by an engine simulation.
Remove those problem data points and keep points in the range [0, Inf].

out = excludedata( Speed, minBSFC, 'Range', [0, Inf] );
f2 = fit( [Speed, Load], minBSFC, 'Lowess', ...
'Normalize', 'on', 'Exclude', out )

Examine the following output:

Locally weighted smoothing linear regression:
       f2(x,y) = lowess (linear) smoothing regression 
       computed from p
       where x is normalized by mean 3443 and std 1187
       and where y is normalized by mean 0.521 and std 0.175
     Coefficients:
       p = coefficient structure

Plot the new fit:

plot( f2, [Speed, Load], minBSFC, 'Exclude', out );
xlabel( 'Speed [RPM]' );
ylabel( 'Load [%]' );
zlabel( 'Minimum BSFC [g/Kwh]' );

Zoom in on the part of the z-axis of interest:
```
set( gca, 'ZLim', [0, max( minBSFC )] );
```
Because you want to operate the engine efficiently, create a contour plot to see the region where the BSFC is low. Use the plot command, and specify the parameter/value pair 'style''Contour'.
```
plot( f2, [Speed, Load], minBSFC, 'Exclude', out,...
 'Style', 'Contour' );
xlabel( 'Speed [RPM]' );
ylabel( 'Load [%]' );
colorbar
```

Create a Table from the Surface

In this exercise, you generate a table from the original data using model f2.

Create variables for the table breakpoints.

speedbreakpoints = linspace( 1000, 5500, 17 );
loadbreakpoints =  linspace( 0.2, 0.8, 13 );

To generate values for the table, evaluate the model over a grid of points.

[tSpeed, tLoad] = meshgrid( speedbreakpoints,...
 loadbreakpoints );
tBSFC = f2( tSpeed, tLoad );

Examine the rows and columns of the table at the command line.
```
tBSFC(1:2:end,1:2:end)
```

Plot the table against the original model. The grid on the model surface shows the table breakpoints.

h = plot( f2 );
set( h, 'EdgeColor', 'none' );
hold on
mesh( tSpeed, tLoad, tBSFC, ...
    'LineStyle', '-', 'LineWidth', 2, 'EdgeColor', 'k', ...
'FaceColor', 'none', 'FaceAlpha', 1 );
hold off
xlabel( 'Speed [RPM]' );
ylabel( 'Load [%]' );
zlabel( 'Minimum BSFC [g/Kwh]' );

Check the table accuracy:

View the difference between the model and the table by plotting the difference between them on a finer grid.
Then, use this difference in prediction accuracy between the table and the model to determine the most efficient table size for your accuracy requirements.

The following code evaluates the model over a finer grid and plots the difference between the model and the table:

[tfSpeed, tfLoad] = meshgrid( linspace( 1000, 5500,...
 8*17+1 ), linspace( 0.2, 0.8, 8*13+1 ) );
tfBSFC_model = f2( tfSpeed, tfLoad );
tfBSFC_table = interp2( tSpeed, tLoad, tBSFC, tfSpeed,...
tfLoad, 'linear' );
tfDiff = tfBSFC_model - tfBSFC_table;

surf( tfSpeed, tfLoad, tfDiff, 'LineStyle', 'none'  );
hold on
mesh( tSpeed, tLoad, zeros( size( tBSFC ) ), ...
    'LineStyle', '-', 'LineWidth', 2, 'EdgeColor', 'k', ...
'FaceColor', 'none', 'FaceAlpha', 1 );
hold off
axis tight
xlabel( 'Speed [RPM]' );
ylabel( 'Load [%]' );
zlabel( 'Difference between model and table [g/Kwh]' );
title( sprintf( 'Max difference: %g', max( abs( tfDiff(:) ) ) ) );

If you have Simulink^® software, you can create a Look Up Table block.

Create a model with a Lookup Table (2-D) block.

simulink
new_system('my_model')
open_system('my_model')
add_block('Simulink/Lookup Tables/Lookup Table (2-D)',...
'my_model/surfaceblock')

Populate the Lookup Table with speed breakpoints, load breakpoints, and a lookup table.

set_param('my_model/surfaceblock',...
'Rowindex','loadbreakpoints',...
'Columnindex','speedbreakpoints',...
'Table','tBSFC');

Examine the populated Lookup Table block.

Biopharmaceutical Drug Interaction Programmatic Surface Fitting Example

Curve Fitting Toolbox provides some example data for an anesthesia drug interaction study. You can use surface fitting methods to fit response surfaces to this data to analyze drug interaction effects.

The following code, using Curve Fitting Toolbox methods, reproduces the interactive surface building with Surface Fitting Tool described in Biopharmaceutical Interactive Surface Fitting Example.

The data is based on the results in the following paper:

Kern SE, Xie G, White JL, Egan TD. Opioid-hypnotic synergy: A response surface analysis of propofol-remifentanil pharmacodynamic interaction in volunteers. Anesthesiology 2004; 100: 1373–81.

Anesthesia is typically at least a two-drug process, consisting of an opioid and a sedative hypnotic. This example uses Propofol and Reminfentanil as drug class prototypes. Their interaction is measured by four different measures of the analgesic and sedative response to the drug combination. Algometry, Tetany, Sedation and Laryingoscopy comprise the four measures of surrogate drug effects at various concentration combinations of Propofol and Reminfentanil.

The following code models the response surfaces for this drug combination:

Load Data

Load the data from file as follows:

data = importdata( 'OpioidHypnoticSynergy.txt' );
Propofol      = data.data(:,1);
Remifentanil  = data.data(:,2);
Algometry     = data.data(:,3);
Tetany        = data.data(:,4);
Sedation      = data.data(:,5);
Laryingoscopy = data.data(:,6);

Create Model Fit Type

Create the model fit type as follows:

ft = fittype( 'Emax*( CA/IC50A + CB/IC50B + alpha*( CA/IC50A ) * ( CB/IC50B ) )^n /(( CA/IC50A + CB/IC50B + alpha*( CA/IC50A ) * ( CB/IC50B ) )^n  + 1 )', ...
    'indep', {'CA', 'CB'}, 'depend', 'z', 'problem', 'Emax' )

Output:

ft = 

     General model:
ft(IC50A,IC50B,alpha,n,Emax,CA,CB) = Emax*...
( CA/IC50A + CB/IC50B + alpha*(CA/IC50A )...
 * ( CB/IC50B ) )^n /(( CA/IC50A + CB/IC50B...
+ alpha*( CA/IC50A ) * ( CB/IC50B ) )^n  + 1 )

Assume Emax = 1

Emax = 1;

Fit a Surface to Algometry

Set fit options as follows:

opts = fitoptions( ft );
opts.Lower = [0 0 -5 -0];
opts.Robust = 'LAR';
opts.StartPoint = [0.00893838724332152 0.706165672266879...
 1 0.746030748284422];

Fit and plot a surface for Algometry:

[f, gof] = fit( [Propofol, Remifentanil], Algometry, ft,...
 opts, 'problem', Emax )
plot( f, [Propofol, Remifentanil], Algometry );

Fit a Surface to Tetany

Fit a surface to Tetany as follows:

[f, gof] = fit( [Propofol, Remifentanil], Tetany, ft, opts, 'problem', Emax )
plot( f, [Propofol, Remifentanil], Tetany );

Fit a Surface to Sedation

Fit a surface to Sedation as follows:

[f, gof] = fit( [Propofol, Remifentanil], Sedation, ft, opts, 'problem', Emax )
plot( f, [Propofol, Remifentanil], Sedation );

Fit a Surface to Laryingoscopy

Fit a surface to Laryingoscopy as follows:

[f, gof] = fit( [Propofol, Remifentanil], Laryingoscopy, ft, opts, 'problem', Emax )
plot( f, [Propofol, Remifentanil], Laryingoscopy );

Generating M-files From Curve Fitting Tool

Curve Fitting Techniques

Recommended Products

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit