I need to make optimization between desired force and displacement curve and simulated one that result from ansys workbench but the RMSE is increase what is the problem?

function error = objectiveFunction(x, desiredForce)
    % Persistent variable to count the number of times the function is called
    persistent ansysCounter;
    if isempty(ansysCounter)
        ansysCounter = 0; % Initialize if empty
    end
    ansysCounter = ansysCounter + 1;
    disp(['ANSYS execution count: ' num2str(ansysCounter)]);

    % Call the external function to handle ANSYS-related processes
    [simulatedForce] = runANSYSProcess(x, ansysCounter);

    % Calculate the error using Root Mean Square Error (RMSE)
    error = sqrt(mean((simulatedForce - desiredForce).^2));  % Scalar value
% Display the calculated RMSE
    disp(['Calculated RMSE: ' num2str(error)]);
end
--------------------------------------------------------------------------------------------------------------
function [simulatedForce] = runANSYSProcess(x, ansysCounter)
    % Extract parameters D1 to D11, PI1 to PI10, H1 to H10 from the input vector x
    D1 = x(1); D2 = x(2); D3 = x(3); D4 = x(4); D5 = x(5);
   D6 = x(6); D7 = x(7); D8 = x(8); D9 = x(9); K99 = x(10); S55 = x(11);
    H1 = x(12); H2 = x(13); H3 = x(14); H4 = x(15); H5 = x(16);
    H6 = x(17); H7 = x(18); H8 = x(19); H9 = x(20); M33 = x(21);
    L40=x(22);

    % Define the path to the base ANSYS journal file
    journalFile ="E:\SPRING\first spring.wbjn";
    if exist(journalFile, 'file') ~= 2
        error('Journal file not found: %s', journalFile);
    end

    % Read the journal content into a string
    journalContent = fileread(journalFile);

    % Replace placeholders in the journal file with parameter values
    journalContent = strrep(journalContent, 'D1', num2str(D1));
    journalContent = strrep(journalContent, 'D2', num2str(D2));
    journalContent = strrep(journalContent, 'D3', num2str(D3));
    journalContent = strrep(journalContent, 'D4', num2str(D4));
    journalContent = strrep(journalContent, 'D5', num2str(D5));
    journalContent = strrep(journalContent, 'D6', num2str(D6));
    journalContent = strrep(journalContent, 'D7', num2str(D7));
    journalContent = strrep(journalContent, 'D8', num2str(D8));
    journalContent = strrep(journalContent, 'D9', num2str(D9));
    journalContent = strrep(journalContent, 'K99', num2str(K99)); 
    journalContent = strrep(journalContent, 'S55', num2str(S55));  
    journalContent = strrep(journalContent, 'H1', num2str(H1));
    journalContent = strrep(journalContent, 'H2', num2str(H2));
    journalContent = strrep(journalContent, 'H3', num2str(H3));
    journalContent = strrep(journalContent, 'H4', num2str(H4));
    journalContent = strrep(journalContent, 'H5', num2str(H5));
    journalContent = strrep(journalContent, 'H6', num2str(H6));
    journalContent = strrep(journalContent, 'H7', num2str(H7));
    journalContent = strrep(journalContent, 'H8', num2str(H8));
    journalContent = strrep(journalContent, 'H9', num2str(H9));
    journalContent = strrep(journalContent, 'M33', num2str(M33)); % Ensure this is being updated
    journalContent = strrep(journalContent, 'L40', num2str(L40)); % Ensure this is being updated

    % Save the modified journal file
    modifiedJournalFile = 'finaljournal.wbjn';
    fid = fopen(modifiedJournalFile, 'w');
    if fid == -1
        error('Failed to open final journal file for writing: %s', modifiedJournalFile);
    end
    fprintf(fid, '%s', journalContent);
    fclose(fid);

    
    % Display the optimized parameter values for debugging
    disp('Optimized Parameters from GA:');
    disp('--------------------------------------------');
    disp(['D1 = ', num2str(D1)]);
    disp(['D2 = ', num2str(D2)]);
    disp(['D3 = ', num2str(D3)]);
    disp(['D4 = ', num2str(D4)]);
    disp(['D5 = ', num2str(D5)]);
    disp(['D6 = ', num2str(D6)]);
    disp(['D7 = ', num2str(D7)]);
    disp(['D8 = ', num2str(D8)]);
    disp(['D9 = ', num2str(D9)]);
    disp(['K99 = ', num2str(K99)]);
    disp(['S55 = ', num2str(S55)]);
    disp(['H1 = ', num2str(H1)]);
    disp(['H2 = ', num2str(H2)]);
    disp(['H3 = ', num2str(H3)]);
    disp(['H4 = ', num2str(H4)]);
    disp(['H5 = ', num2str(H5)]);
    disp(['H6 = ', num2str(H6)]);
    disp(['H7 = ', num2str(H7)]);
    disp(['H8 = ', num2str(H8)]);
    disp(['H9 = ', num2str(H9)]);
    disp(['M33 = ', num2str(M33)]);
     disp(['L40 = ', num2str(L40)]);

    disp('--------------------------------------------');
    % Run ANSYS using the modified journal file
    status = system('"C:\Program Files\ANSYS Inc\v242\Framework\bin\Win64\runwb2.exe" -x -R "finaljournal.wbjn"');
    if status ~= 0
        error('ANSYS execution failed. Status code: %d', status);
    end

    % Wait for ANSYS to finish by checking if the output file exists
    csvFile ="E:\SPRING\first spring.csv"; % Path to the ANSYS output CSV file
    if ~isfile(csvFile)
        error('ANSYS output file not found: %s', csvFile);
    end

    % Load the force and displacement data from the CSV file generated by ANSYS
    simulatedForce = readmatrix(csvFile, 'Range', 'BA8:CC8');
    if isempty(simulatedForce)
        error('Failed to load force data from file: %s', csvFile);
    end

desiredDeformation = [ 0.348347797, 1.173153205, 2.21822933, 3.812862584, 4.803884609, 5.90554205, 6.953464755, 7.89172824, 8.664936558, 9.714926366, 10.60098324, 11.59845866, 12.48718235, 13.32165075, 14.26752652, 15.15919716, 15.88505542, 16.66788445, 17.34050795, 18.01377124, 18.63121851, 19.13718778, 19.75680325, 20.43350892, 20.88563395, 21.22577958, 21.79170524, 22.13235055, 22.41673682];


    % Plot optimized force-deformation curve every 100 ANSYS executions
    if mod(ansysCounter, 1) == 0
        figure(1); % Open or create figure window 1
        hold on; % Retain the current plot to overlay new data
        plot(desiredDeformation, simulatedForce, '-o', ...
             'DisplayName', ['ANSYS Run ' num2str(ansysCounter)], 'LineWidth', 1);
    end
end
------------------------------------------------------------------------------------------------------------------------------------------------------------
% Main script for running the optimization loop
clear objectiveFunction; % Clear persistent variables in the objective function

% Define desired deformation and force data as row vectors
desiredDeformation = [ 0.348347797, 1.173153205, 2.21822933, 3.812862584, 4.803884609, 5.90554205, 6.953464755, 7.89172824, 8.664936558, 9.714926366, 10.60098324, 11.59845866, 12.48718235, 13.32165075, 14.26752652, 15.15919716, 15.88505542, 16.66788445, 17.34050795, 18.01377124, 18.63121851, 19.13718778, 19.75680325, 20.43350892, 20.88563395, 21.22577958, 21.79170524, 22.13235055, 22.41673682];
desiredForce =[-188.7501842, -200.5384007, -217.2529154, -232.9123719, -249.6333495, -267.3318313, -289.9889256, -309.6876558, -325.4440548, -344.1394294, -368.7967721, -392.4507591, -423.0506813, -454.6474963, -480.288806, -505.9365786, -530.6133097, -558.2548678, -588.8806413, -617.525555, -639.2439221, -658.994355, -687.6457315, -717.2810752, -740.9996906, -765.7216616, -792.3986411, -815.1397523, -836.8968964];
% Plot desired force-deformation curve
figure(1); % Create or select figure 1
clf; % Clear current figure
hold on; % Retain plots for overlaying
plot(desiredDeformation, desiredForce, '-x', 'DisplayName', 'Desired Deformation', 'LineWidth', 1.5);
xlabel('Deformation');
ylabel('Force');
title('Comparison between Optimized and Desired Force-Deformation Curves');
legend('show');
grid on;

% Define number of variables and bounds
nVars = 22;  % Number of optimization variables

% Lower and upper bounds for the optimization variables
lb =[15, 20, 25, 30, 40, 45, 40, 30, 25, 20, 15, 7,7,7,7,7,7,7,7,7,7,7];
ub = [45, 50, 55, 60, 75, 82, 75, 65, 60, 55, 48, 15,15,15,15,15,15,15,15,15,15];
initialCondition = [18, 25, 30, 35, 45, 53.9, 45, 35, 30, 25, 18, 8,8.5,9,9.5,10,9.5,9,8.5,8,7,7];


options = optimoptions('ga', ...
    'PopulationSize', 200, ...  % Size of the population
    'CrossoverFraction', 0.8, ...  % Fraction of the population to cross-over
    'MutationFcn', {@mutationadaptfeasible}, ...  % Mutation function
    'EliteCount', 20, ...  % Number of elite solutions to preserve
    'SelectionFcn', @selectiontournament, ...  % Tournament selection
    'MaxGenerations', Inf, ...  % Maximum number of generations
    'MaxStallGenerations', Inf, ...  % Maximum number of generations with no improvement
    'OutputFcn', @gaOutputFcn, ...  % Custom output function
    'Display', 'iter', ...  % Display output at each iteration
    'PlotFcn', {@gaplotbestf, @gaplotstopping}, ... % Plot best function and stopping criteria
    'InitialPopulation', initialCondition);  % Provide the initial condition as one individual


% Run the Genetic Algorithm
[x_opt, fval] = ga(@(x) objectiveFunction(x, desiredForce), nVars, [], [], [], [], lb, ub, @nonlcon, options);

% Display the optimized parameters
disp('Optimized Parameters:');
for i = 1:11
    fprintf('D%d: %.2f\n', i, x_opt(i));
end
for i = 12:20
    fprintf('H%d: %.2f\n', i, x_opt(i));
end
for i = 21
    fprintf('M%d: %.2f\n', i, x_opt(i));
end
for i = 22
    fprintf('L%d: %.2f\n', i, x_opt(i));
end

% Load the final optimized data from ANSYS
optimizedForce = readmatrix("E:\SPRING\first spring.csv", 'Range', 'BA8:CC8'); % Updated range

% Plot the final optimized force-deformation curve
figure(1);
plot(desiredDeformation, optimizedForce, '-o', 'DisplayName', 'Final Optimized', 'LineWidth', 1.5);

% Update the legend and display the plot
legend('show');
drawnow;
hold off;