Exponenta Banner

Тематическое исследование бензина

Этот пример использует:

Скрипт Open Live Script

В этом примере показано, как автоматически сгенерировать mbcmodel проект для тематического исследования бензина с помощью функций командной строки. Во-первых, вы создаете и загружаете данные в проект. Затем вы создаете план тестирования, граничные модели и ответы. Затем вы удаляете выбросы данных и создаете модели функции. Наконец, вы делаете модель 2D этапа и строите крутящий момент по сравнению с искрой. Пример использует DIVCP_Main_DoE_Data.mat файл в mbctraining папка.

Создайте новый mbcmodel Проект

Используйте CreateProject функция, чтобы создать проект.

project = mbcmodel.CreateProject;

Загрузите данные в проект

Используйте функции, чтобы сгруппировать данные в тесты, удалить неправильные данные и удалить тесты, которые не имеют достаточного количества данных, чтобы подбирать локальные модели.

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 );

Создайте план тестирования

Задайте входные параметры для плана тестирования.

LocalInputs = mbcmodel.modelinput('Symbol','S',...
    'Name','SPARK',...
    'Range',[0 50]);
GlobalInputs = mbcmodel.modelinput('Symbol',{'N','L','ICP','ECP'},...
    'Name',{'SPEED','LOAD','INT_ADV','EXH_RET'},...
    '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 );

Создайте граничные модели

Создайте глобальную граничную модель и добавьте его в дерево плана тестирования. 

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);

% 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.
testplan.Boundary.Global
ans = 
  Tree with properties:

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

Создайте ответы

Создайте модели ответа для крутящего момента, исчерпайте температурную и остаточную часть. Используйте локальную полиномиальную модель сплайна для крутящего момента. Для выхлопной температурной и остаточной части используйте локальную полиномиальную модель с данной величиной.

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

Удалите локальные выбросы

Удалите данные если abs(studentized residuals) > 3. Gasoline_project использует различный процесс, чтобы решить который выбросы удалить.

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);
    end
    % get local model for test and look at summary statistics
    mdl = ModelForTest(LocalBTQ,tn);
    if ~strcmp(mdl.Status,'Not fitted')
        LocalStats = SummaryStatistics(mdl);
    end
end

Обновите функции ответа.

UpdateResponseFeatures(LocalBTQ);

Удалите Точки Где MBT <0 или MBT> 60

Удалите точки, где максимальный момент привода (MBT) меньше 0 или больше, чем 60.

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

Создайте альтернативные модели функции ответа

Составьте список этих альтернативных моделей функции ответа. Выберите лучшую модель, на основе информационных Критериев Akaike (AICc).

  • Квадратичный

  • Кубический

  • RBF с областью значений центров

  • Полином-RBF с областью значений центров

Получите базовую модель. Вы используете это, чтобы создать другие модели.

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

Сделайте квадратичную модель, которая использует Minimize PRESS соответствовать. Добавьте его в список.

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

Сделайте кубическую модель и добавьте его в список.

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

Сделайте модели RBF с областью значений центров. Максимальный номер центров определяется в центральном алгоритме выбора.

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

Сделайте модели Polynomial-RBF с областью значений центров.

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;
end

Заставьте альтернативные модели для каждого ответа показать и выбрать лучшую модель, на основе AICc.

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

Измените модели функции ответа

Получите альтернативные ответы для модели узла. Измените модели с помощью ступенчатой регрессии.

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

Получите пошаговую статистику.

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

Используйте PRESS найти, что самое благоприятное условие изменяется. Переключите пошаговое состояние термина с индексом.

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

Получите статистическую величину фактора инфляции отклонения (VIF).

VIF = MultipleVIF(knot_model)
VIF = 11×1

    3.1290
    1.2918
    1.6841
    1.1832
    1.3230
    2.6617
    1.6603
    1.3306
    1.2856
    1.4317
      ⋮

Получите RMSE.

RMSE = SummaryStatistics(knot_model, 'RMSE')
RMSE = 5.1578

Измените модель в Полином-RBF, имеющий до 10 центров.

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.
[new_knot_model,S] = fit(new_knot_model);

Если нет никаких проблем с изменениями, не обновляют ответ. В противном случае продолжите использовать исходную модель.

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

Сделайте модель 2D этапа для крутящего момента

Сделайте модель 2D этапа для ответа крутящего момента.

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

    0.8319    4.6117

Постройте модель 2D этапа крутящего момента против Spark

Постройте крутящий момент по сравнению с искрой для теста 5.

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

Создайте другие ответы

Выполните эти шаги, чтобы создать модели ответа для выхлопной температурной и остаточной части.

  • Скопируйте выбросы с модели крутящего момента.

  • Сделайте альтернативные модели для каждой функции ответа.

  • Сделайте модель 2D этапа без оценки наибольшего правдоподобия (MLE).

Исчерпайте температурный ответ

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

   10.5648   27.9941

Остаточный дробный ответ

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

    0.0824    0.5596


if isgraphics(fig)
    % delete figure made during example
    delete(fig)
end

Смотрите также

| | | | | |

Документация Model-Based Calibration Toolbox

Поддержка