Этот пример показывает, как автоматически сгенерировать mbcmodel проект для тематического исследования бензина с помощью инструментов командной строки в Model-Based Calibration Toolbox™.
Требует DIVCP_Main_DoE_Data.xls от mbctraining папки.
mbcmodel
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 );
Создайте глобальную граничную модель. CreateBoundary не добавляет граничную модель в дерево.
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. testplan.Boundary.Global
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]
Создайте модели ответа для крутящего момента, исчерпайте температурную и остаточную часть * Использование локальная полиномиальная модель сплайна для крутящего момента. * Используйте локальную полиномиальную модель с данной величиной для выхлопной температурной и остаточной части
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 невязки)> 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);
Starting parallel pool (parpool) using the 'local' profile ... Connected to the parallel pool (number of workers: 12).
knot = LocalBTQ.ResponseFeatures(1); PointsToRemove = knot.DoubleResponseData<0 | knot.DoubleResponseData>60; knot.RemoveOutliers(PointsToRemove);
Составьте список моделей кандидата и выберите лучшее на основе AICc.
Квадратичный
Кубический
RBF с областью значений центров
Полином-RBF с областью значений центров
Получите базовую модель. Вы используете это, чтобы создать другие модели.
rf = LocalBTQ.ResponseFeatures(1); BaseModel = rf(1).Model;
Сделайте квадратичную модель, которой использование Минимизирует НАЖАТИЕ, чтобы соответствовать и добавить его в список.
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) m.FitAlgorithm.WidthAlgorithm.NestedFitAlgorithm.CenterSelectionAlg.MaxCenters = Centers(i); 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);
Используйте НАЖАТИЕ, чтобы найти, что самое благоприятное условие изменяет, и переключает пошаговое состояние термина с этим индексом.
[bestPRESS, ind] = min(stepwise_stats(:,4)); [stepwise_stats,knot_model] = StepwiseRegression(knot_model, ind);
Получите статистическую величину VIF.
VIF = MultipleVIF(knot_model)
VIF = 3.1290 1.2918 1.6841 1.1832 1.3230 2.6617 1.6603 1.3306 1.2856 1.4317 2.6550
Получите 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
doMLE = true; MakeHierarchicalResponse( LocalBTQ, doMLE ); % Look at the Local and Two-Stage RMSE. BTQ_RMSE = SummaryStatistics( LocalBTQ, {'Local RMSE', 'Two-Stage RMSE'} )
BTQ_RMSE = 0.8319 4.6118
Постройте 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
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 = 10.5648 27.9916
Ответ RESIDFRAC
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'} ) if isgraphics(fig) % delete figure made during example delete(fig) end
RESIDFRAC_RMSE = 0.0824 0.5596