we have input data fuel and output shaft speed of turbojet engine , and extract model using % Load data

%load('data.mat'); % Contains input (fuel) and output (shaft speed)

load('C:\Users\HP\Downloads\Data.mat');

X = Output2; % System output (shaft speed)

U = input; % System input (fuel)

% Normalize data

X_normalized = (X - mean(X)) / std(X);

U_normalized = (U - mean(U)) / std(U);

% Prepare time-series data for NARX

inputSeries = tonndata(U_normalized, false, false);

targetSeries = tonndata(X_normalized, false, false);

% Define NARX parameters

inputDelays = 1:2; % Past inputs (t-1, t-2)

feedbackDelays = 1:2; % Past outputs (t-1, t-2)

hiddenLayerSize = 10; % Number of neurons in the hidden layer

% Create NARX network

narxNet = narxnet(inputDelays, feedbackDelays, hiddenLayerSize);

% Prepare input and target data

[Xs, Xi, Ai, Ts] = preparets(narxNet, inputSeries, {}, targetSeries);

% Divide data into training, validation, and test sets

narxNet.divideFcn = 'divideblock';

narxNet.divideParam.trainRatio = 70/100;

narxNet.divideParam.valRatio = 15/100;

narxNet.divideParam.testRatio = 15/100;

% Train NARX network

narxNet.trainFcn = 'trainlm'; % Levenberg-Marquardt training

narxNet.trainParam.epochs = 100; % Maximum number of epochs

narxNet.trainParam.min_grad = 1e-10; % Minimum gradient

narxNet.trainParam.max_fail = 10; % Maximum validation failures

% Train the network

[narxNet, tr] = train(narxNet, Xs, Ts, Xi, Ai);

% Simulate the network

Y = sim(narxNet, Xs, Xi, Ai);

% Convert output back to matrix

Y = cell2mat(Y);

Ts = cell2mat(Ts); % Convert target series to matrix for comparison

% Compute errors

RMSE_narx = sqrt(mean((Ts - Y).^2)); % Fixed syntax error

MAE_narx = mean(abs(Ts - Y));

corr_coeff_narx = corrcoef(Ts, Y);

% Display results

disp(['NARX RMSE: ', num2str(RMSE_narx)]);

disp(['NARX MAE: ', num2str(MAE_narx)]);

disp(['NARX Correlation Coefficient: ', num2str(corr_coeff_narx(1,2))]);

% Plot results

figure;

subplot(2,1,1);

plot(Ts, 'b', 'LineWidth', 2); hold on;

plot(Y, 'r--', 'LineWidth', 2);

legend('Actual Data', 'NARX Predicted');

xlabel('Time Step');

ylabel('Shaft Speed');

title('Turbojet Engine: NARX Model Validation');

grid on;

% Residual Analysis

subplot(2,1,2);

plot(Ts - Y, 'k', 'LineWidth', 1.5);

xlabel('Time Step');

ylabel('Residual Error');

title('Residual Analysis');

grid on;

% 2 Fit Linear Regression Model

linModel = fitlm(U, X);

% Predict Shaft Speed

Y_lin = predict(linModel, U);

% Compute Errors

RMSE_lin = sqrt(mean((X - Y_lin).^2));

MAE_lin = mean(abs(X - Y_lin));

corr_lin = corrcoef(X, Y_lin);

% Display Results

disp(['Linear Regression RMSE: ', num2str(RMSE_lin)]);

disp(['Linear Regression MAE: ', num2str(MAE_lin)]);

disp(['Linear Regression Correlation: ', num2str(corr_lin(1,2))]);

% Define ARX Model Orders

na = 2; % Past output orders

nb = 2; % Past input orders

nk = 1; % Input delay

% Create IDDATA object

data = iddata(X, U, 1); % Assumes Ts = 1 (modify if needed)

% Identify ARX Model

sys_arx = arx(data, [na nb nk]);

% Predict Output

Y_arx = predict(sys_arx, data);

% Convert Predicted Output to Array

Y_arx = Y_arx.OutputData;

% Compute Errors

RMSE_arx = sqrt(mean((X - Y_arx).^2));

MAE_arx = mean(abs(X - Y_arx));

corr_arx = corrcoef(X, Y_arx);

% Display Results

disp(['ARX RMSE: ', num2str(RMSE_arx)]);

disp(['ARX MAE: ', num2str(MAE_arx)]);

disp(['ARX Correlation: ', num2str(corr_arx(1,2))]);

% Plot Results

figure;

subplot(2,1,1);

plot(X, 'b', 'LineWidth', 2); hold on;

plot(Y_arx, 'r--', 'LineWidth', 2);

legend('Actual Data', 'ARX Predicted');

xlabel('Time Step');

ylabel('Shaft Speed');

title('Turbojet Engine: ARX Model Validation');

grid on;

% Residual Analysis

subplot(2,1,2);

plot(X - Y_arx, 'k', 'LineWidth', 1.5);

xlabel('Time Step');

ylabel('Residual Error');

title('Residual Analysis');

grid on;

Guide how to make MPC