This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English verison of the page.

Note: This page has been translated by MathWorks. Please click here
To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

Gasoline Case Study

This example shows how to automatically generate an mbcmodel project for the gasoline case study using the command-line tools in Model-Based Calibration Toolbox™ .

Requires DIVCP_Main_DoE_Data.xls from mbctraining folder.

Create a New mbcmodel Project

project = mbcmodel.CreateProject;

Load Data into Project

  • Group data into tests and add some filters.

  • Add filters to remove bad data.

  • Remove tests which do not have sufficient data to fit local models.

datafile = fullfile( mbcpath, 'mbctraining', 'DIVCP_Main_DoE_Data.mat' );
data = CreateData( project, datafile );

BeginEdit( data );

% Group data by test number GDOECT.
DefineTestGroups( data, 'GDOECT', 0.5, 'GDOECT', false );
% Get rid of data which is probably unstable.
AddFilter( data, 'RESIDFRAC < 35' );
AddFilter( data, 'AFR > 14.25' );
% Get rid of the tests that are too small.
AddTestFilter( data, 'length(BTQ) > 4' );

% Commit the changes to the project.
CommitEdit( data );

Build Test Plan

Define Inputs for test plan.

LocalInputs = mbcmodel.modelinput('Symbol','S',...
    'Range',[0 50]);
GlobalInputs = mbcmodel.modelinput('Symbol',{'N','L','ICP','ECP'},...
    'Range',{[500 6000],[0.0679    0.9502],[-5 50],[-5 50]});
% Create test plan.
testplan = CreateTestplan( project, {LocalInputs,GlobalInputs} );
% Attach data to the test plan.
AttachData( testplan, data );

Build Boundary Models

Create a global boundary model. CreateBoundary does not add the boundary model to the tree.

B = CreateBoundary(testplan.Boundary.Global,'Star-shaped');
% Add boundary model to the test plan. The boundary model is fitted when it
% is added to the boundary model tree. The boundary model is included in
% the best boundary model for the tree by default.
% All inputs are used in the boundary model by default.
B = Add(testplan.Boundary.Global,B);

% Now make a boundary model using only speed and load and add to the
% boundary tree.
B.ActiveInputs = [true true false false];
B = Add(testplan.Boundary.Global,B);
% Look at the global boundary tree.
ans = 

  Tree with properties:

         Data: [189x4 double]
       Models: {[1x1 mbcboundary.Model]  [1x1 mbcboundary.Model]}
    BestModel: [1x1 mbcboundary.Boolean]
       InBest: [1 1]
     TestPlan: [1x1 mbcmodel.testplan]

Build Responses

Build response models for torque, exhaust temperature and residual fraction * Use a local polynomial spline model for torque. * Use a local polynomial model with datum for exhaust temperature and residual fraction

LocalTorque  = mbcmodel.CreateModel('Local Polynomial Spline',1);
LocalTorque.Properties.LowOrder = 2;
% Use the default global model.
GlobalModel = testplan.DefaultModels{2};
% make exhaust temperature and residual fraction models
LocalPoly  = mbcmodel.CreateModel('Local Polynomial with Datum',1);

Remove Local Outliers

Remove data if abs(studentized residuals) > 3. Note that a different process was used in the project Gasoline_project to decide which outliers to remove.

TQ_response = testplan.Responses(1);
numTests = TQ_response.NumberOfTests;
LocalBTQ = TQ_response.LocalResponses;
for tn = 1:numTests
    % Find observations with studentized residuals greater than 3
    studentRes = DiagnosticStatistics( LocalBTQ, tn, 'Studentized residuals' );
    potentialOut  = abs( studentRes )> 3;
    if any(potentialOut)
        % Don't update response feature models until end of loop
        RemoveOutliersForTest( LocalBTQ, tn, potentialOut , false);
    % get local model for test and look at summary statistics
    mdl = ModelForTest(LocalBTQ,tn);
    if ~strcmp(mdl.Status,'Not fitted')
        LocalStats = SummaryStatistics(mdl);

Update response features.

Starting parallel pool (parpool) using the 'local' profile ...
connected to 12 workers.

Remove Points Where MBT<0 or MBT>60

knot = LocalBTQ.ResponseFeature(1);
PointsToRemove = knot.DoubleResponseData<0 | knot.DoubleResponseData>60;

Create Alternative Response Feature Models

Make a list of candidate models and select the best based on AICc.

  • Quadratic

  • Cubic

  • RBF with a range of centers

  • Polynomial-RBF with a range of centers

Get the base model. You use this to create the other models.

rf = LocalBTQ.ResponseFeature(1);
BaseModel = rf(1).Model;

Make a quadratic model that uses Minimize PRESS to fit, and add it to the list.

m = BaseModel.CreateModel('Polynomial');
m.Properties.Order = [2 2 2 2];
m.FitAlgorithm = 'Minimize PRESS';
mlist = {m};

Make a cubic model and add it to the list.

m.Properties.Order = [3 3 3 3];
m.Properties.InteractionOrder = 2;
mlist{2} = m;

Make RBF models with a range of centers. The maximum number of centers is set in the center selection algorithm.

m = BaseModel.CreateModel('RBF');
Centers = [50 80];
Start = length(mlist);
mlist = [mlist cell(size(Centers))];
for i = 1:length(Centers)
    m.FitAlgorithm.WidthAlgorithm.NestedFitAlgorithm.CenterSelectionAlg.MaxCenters = Centers(i);
    mlist{Start+i} = m;

Make Polynomial-RBF models with a range of centers.

m = BaseModel.CreateModel('Polynomial-RBF');
m.Properties.Order = [2 2 2 2];
Start = length(mlist);
mlist = [mlist cell(size(Centers))];
for i = 1:length(Centers)
    % Maximum number of centers is set in the nested fit algorithm
    m.FitAlgorithm.WidthAlgorithm.NestedFitAlgorithm.MaxCenters = Centers(i);
    mlist{Start+i} = m;

Make alternative models for each response feature and select the best model based on AICc.

criteria = 'AICc';
CreateAlternativeModels( LocalBTQ, mlist, criteria );

Alter Response Feature Models

Get the alternative responses for knot and alter models using stepwise regression.

knot = LocalBTQ.ResponseFeature(1);
AltResponse = knot.AlternativeResponses(1);

Get the stepwise statistics.

knot_model = AltResponse.Model;
[stepwise_stats,knot_model] = StepwiseRegression(knot_model);

Use PRESS to find the best term to change, and toggle the stepwise status of the term with this index.

[bestPRESS, ind] = min(stepwise_stats(:,4));
[stepwise_stats,knot_model] = StepwiseRegression(knot_model, ind);

Get a VIF statistic.

VIF = MultipleVIF(knot_model)


Get the RMSE.

RMSE = SummaryStatistics(knot_model, 'RMSE')


Change the model to a Polynomial-RBF with a maximum of 10 centers.

new_knot_model = knot_model.CreateModel('Polynomial-RBF');
new_knot_model.Properties.Order = [1 1 1 1];
new_knot_model.FitAlgorithm.WidthAlgorithm.NestedFitAlgorithm.MaxCenters = 10;
% Fit the model with current data.
[S,new_knot_model] = new_knot_model.Fit;

If there were no problems with the changes then update the response, otherwise you will continue to use the original model.

if strcmp(new_knot_model.Status,'Fitted')
    new_RMSE = SummaryStatistics(new_knot_model,'RMSE')
    % Update the response with the new model.
new_RMSE =


Make the Two-stage Model for Torque

doMLE = true;
MakeHierarchicalResponse( LocalBTQ, doMLE );
% Look at the Local and Two-Stage RMSE.
BTQ_RMSE = SummaryStatistics( LocalBTQ, {'Local RMSE', 'Two-Stage RMSE'} )

    0.8319    4.6117

Plot the Two-stage Model of Torque Against Spark

Plot the 5th test

testToPlot = 5;
BTQInputData = TQ_response.DoubleInputData(testToPlot);
BTQResponseData = TQ_response.DoubleResponseData(testToPlot);
BTQPredictedValue = TQ_response.PredictedValue( BTQInputData );
fig = figure;
plot( BTQInputData(:,1), BTQResponseData, 'o', BTQInputData(:,1), BTQPredictedValue, '-' );
xlabel( 'spark' );
ylabel( 'torque' );
title( 'Test 5' );
grid on

Build Other Responses

Build models for exhaust temperature and residual fraction.

  • Copy outliers from torque model

  • Make alternative models for each response feature

  • Make two-stage model without MLE

EXTEMP Response

EXTEMP = testplan.Responses(2).LocalResponses(1);
CreateAlternativeModels( EXTEMP,mlist, criteria );
MakeHierarchicalResponse( EXTEMP, false );
EXTEMP_RMSE = SummaryStatistics( EXTEMP, {'Local RMSE', 'Two-Stage RMSE'} )

   10.5648   27.9919


RESIDFRAC = testplan.Responses(3).LocalResponses(1);
CreateAlternativeModels( RESIDFRAC,mlist, criteria );
ok = MakeHierarchicalResponse( RESIDFRAC, false );
RESIDFRAC_RMSE = SummaryStatistics( RESIDFRAC, {'Local RMSE', 'Two-Stage RMSE'} )

if isgraphics(fig)
    % delete figure made during example

    0.0824    0.5595

Was this topic helpful?