Error while running my objective function by fmincon with nonlinear constraint for calculate optimization coordination overcurrent relay

i have writing code with AI and this happened:

Error using ^ (line 52)

Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To operate on each element of the matrix

individually, use POWER (.^) for elementwise power.

Error in

objective_fun (line 15)

denominator = I_norm^b - 1;

^

Error in

main_optimization>@(x)objective_fun(x,R_pairs,numRelays,a,b,c,V_i,I_f_i,PS_max_i) (line 71)

objectiveFunction = @(x) objective_fun(x, R_pairs, numRelays, a, b, c, V_i, I_f_i, PS_max_i);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in

objfunEvaluator (line 5)

fval = feval(Objfun, x, self.FunArgs.AdditionalParameters{:});

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in

OptimFunctions/objectiveFirstEval (line 655)

[fval, grad, hess] = self.ObjectiveFunAndGrad(self,self.FunFcn{3},...

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in

fmincon (line 505)

[initVals.f,initVals.g,HESSIAN,funObj] = funObj.objectiveFirstEval( ...

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in

main_optimization (line 75)

[x_optimized, fval, exitflag, output] = fmincon(objectiveFunction, x0, A, b, Aeq, beq, lb, ub, nonlinearConstraints, options);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue.

29 obj.task();

Operation terminated by user during matlab.internal.doc.ui.java.initializeHelpSystem (line 3)

In

matlab.internal.doc.ui.setupForHelpUI (line 10)

matlab.internal.doc.ui.java.initializeHelpSystem;

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In

matlab.internal.doc.ui.CshDocPageHandler/openBrowserForDocPage (line 15)

matlab.internal.doc.ui.setupForHelpUI;

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In

matlab.internal.doc.ui.DocPageHandler/openBrowser (line 4)

success = openBrowserForDocPage(obj, url);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In

matlab.internal.doc.ui.DocPageLauncher/openDocPage (line 31)

success = openBrowser(obj.Handler, getNavigationUrl(activePage));

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In

doc (line 76)

launcher.openDocPage;

^^^^^^^^^^^^^^^^^^^^

In

helpPopup (line 18)

doc(topic);

^^^^^^^^^^

In

matlab.lang.internal.introspective.errorDocCallback>popTopicHelp (line 54)

helpPopup(topic);

^^^^^^^^^^^^^^^^

In

matlab.lang.internal.introspective.errorDocCallback (line 37)

if hasFileMarker && popTopicHelp(topicWithoutLocalFunction)

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

i have tried matlab ai but it didn't helped, my problem is still unsolved

this is the code:

% main_optimization.m
% The main script for solving protective relay coordination optimization problems.
% MATLAB version: R2025a
clear; clc; close all;

%% 1. System Parameter Definition
% Suppose we have 5 relays, where R_p and R_c are defined as
% the primary and backup relay pairs (i,j).
numRelays = 5;
numScenarios = 3; % Number of fault scenarios

% Define the primary-backup (S) relationship in matrix or cell array form
% Example: (primary, backup) pairs for scenarios 1, 2, and 3
R_pairs = cell(numScenarios, 1);
R_pairs{1} = [1, 2; 3, 4]; % Scenario 1: (Relay 1, Relay 2), (Relay 3, Relay 4)
R_pairs{2} = [1, 5; 4, 2]; % Scenario 2: (Relay 1, Relay 5), (Relay 4, Relay 2)
R_pairs{3} = [3, 5];       % Scenario 3: (Relay 3, Relay 5)

% Input known parameters for each relay
% V_i, I_f_i, and PS_max_i for each relay (numRelays)
V_i = [10; 12; 15; 11; 13];
I_f_i = [500; 450; 600; 550; 480];
PS_max_i = [5; 6; 7; 5.5; 6.5];

% Constants for the problem
a = 0.02;
b = 0.14;
c = 0;
CTI_min = 0.2;
t_min = 0.02;
TMS_min = 0.01;
TMS_max = 1;

%% 2. Defining Optimization Variables and Constraints
% Optimization variables: x = [TMS_1, ..., TMS_N, lambda_1, ..., lambda_N, xi_1, ..., xi_N]
% N = numRelays
numVarPerRelay = 3; % TMS, lambda, xi

% Lower and Upper Limits (lb, ub)
% The arrangement of lb and ub must be the same as the arrangement of x
lb_TMS = zeros(numRelays, 1) + TMS_min; % or directly lb_TMS = TMS_min * ones(numRelays, 1);
ub_TMS = zeros(numRelays, 1) + TMS_max;
lb_lambda = zeros(numRelays, 1);
ub_lambda = PS_max_i / 0.9;
lb_xi = zeros(numRelays, 1);
ub_xi = PS_max_i;

lb = [lb_TMS; lb_lambda; lb_xi];
ub = [ub_TMS; ub_lambda; ub_xi];

x0 = (lb + ub) / 2; % Initial guess in the middle of the boundary

% Since there are no linear constraints, set A, b, Aeq, beq to zero.
A = [];
b = [];
Aeq = [];
beq = [];

%%3. Calling the fmincon solver
% We use optimoptions to select the 'interior-point' algorithm and
% set other preferences
options = optimoptions('fmincon', 'Algorithm', 'interior-point', ...
                       'Display', 'iter-detailed', ...
                       'UseParallel', false); % Utilize CPU cores if available

% Objective function and constraint definition
% Note that we use the anonymous function (@) to pass additional parameters
% (such as a, b, c, CTI_min, etc.) to the objective function and constraints.
objectiveFunction = @(x) objective_fun(x, R_pairs, numRelays, a, b, c, V_i, I_f_i, PS_max_i);
nonlinearConstraints = @(x) nonlin_con(x, R_pairs, CTI_min, t_min, numRelays, a, b, c, V_i, I_f_i, PS_max_i);

% Calling fmincon
[x_optimized, fval, exitflag, output] = fmincon(objectiveFunction, x0, A, b, Aeq, beq, lb, ub, nonlinearConstraints, options);

%% 4. Display the result
if exitflag > 0
    fprintf('\nOptimization is successfully finished!\n');
    
    % Extract the result
    optimized_TMS = x_optimized(1:numRelays);
    optimized_lambda = x_optimized(numRelays+1 : 2*numRelays);
    optimized_xi = x_optimized(2*numRelays+1 : 3*numRelays);
    
    fprintf('The value of TMS optimized: \n');
    disp(optimized_TMS);
    fprintf('The value of lambda optimized: \n');
    disp(optimized_lambda);
    fprintf('The value of xi optimized: \n');
    disp(optimized_xi);
    fprintf('Value of minimum objective function (total operation time): %.4f\n', fval);
else
    fprintf('\nOptimization unsuccessfully finished. Exit Flag: %d\n', exitflag);
    disp(output);
end

% objective_fun.m
% Objective function for fmincon.

function total_time = objective_fun(x, R_pairs, numRelays, a, b, c, V_i, I_f_i, PS_max_i)
    % Separate optimized variabel from vector x
    TMS = x(1:numRelays);
    lambda = x(numRelays+1 : 2*numRelays);
    xi = x(2*numRelays+1 : 3*numRelays);

    % Calculate PS and t for every relay
    PS = min(lambda * V_i + xi, PS_max_i);
    I_norm = I_f_i / PS;
    
    % Avoid zero dividen or negative value inside algorithm
    denominator = (I_norm)^b - 1;
    t_i = TMS * (a / denominator) + c;
    % Change NaN or Inf if any
    t_i(denominator <= 0) = inf;

    % Calculate total operation time from relay's pair
    total_time = 0;
    for s = 1:length(R_pairs)
        scenario_pairs = R_pairs{s};
        for k = 1:size(scenario_pairs, 1)
            i_p = scenario_pairs(k, 1); % Primary relay
            i_c = scenario_pairs(k, 2); % Backup Relay
            
            % Pastikan indeks valid
            if i_p > numRelays || i_c > numRelays
                continue; % Pass if the index is invalid
            end

            total_time = total_time + t_i(i_p) + t_i(i_c);
        end
    end
end

Answer 1

[Matt J]

You have a parameter b defined here

b = 0.14;

but then you overwrite it here:

A = [];
b = [];
Aeq = [];
beq = [];

Also, review the operators in your code to make sure you are not using P/Q when you mean P./Q or P^Q when you mean P.^Q.

Answer 2

[Torsten]

So far, it seems to work. The constraint function is missing.

% main_optimization.m
% The main script for solving protective relay coordination optimization problems.
% MATLAB version: R2025a
clear; clc; close all;
%% 1. System Parameter Definition
% Suppose we have 5 relays, where R_p and R_c are defined as
% the primary and backup relay pairs (i,j).
numRelays = 5;
numScenarios = 3; % Number of fault scenarios
% Define the primary-backup (S) relationship in matrix or cell array form
% Example: (primary, backup) pairs for scenarios 1, 2, and 3
R_pairs = cell(numScenarios, 1);
R_pairs{1} = [1, 2; 3, 4]; % Scenario 1: (Relay 1, Relay 2), (Relay 3, Relay 4)
R_pairs{2} = [1, 5; 4, 2]; % Scenario 2: (Relay 1, Relay 5), (Relay 4, Relay 2)
R_pairs{3} = [3, 5];       % Scenario 3: (Relay 3, Relay 5)
% Input known parameters for each relay
% V_i, I_f_i, and PS_max_i for each relay (numRelays)
V_i = [10; 12; 15; 11; 13];
I_f_i = [500; 450; 600; 550; 480];
PS_max_i = [5; 6; 7; 5.5; 6.5];
% Constants for the problem
a = 0.02;
b = 0.14;
c = 0;
CTI_min = 0.2;
t_min = 0.02;
TMS_min = 0.01;
TMS_max = 1;
%% 2. Defining Optimization Variables and Constraints
% Optimization variables: x = [TMS_1, ..., TMS_N, lambda_1, ..., lambda_N, xi_1, ..., xi_N]
% N = numRelays
numVarPerRelay = 3; % TMS, lambda, xi
% Lower and Upper Limits (lb, ub)
% The arrangement of lb and ub must be the same as the arrangement of x
lb_TMS = zeros(numRelays, 1) + TMS_min; % or directly lb_TMS = TMS_min * ones(numRelays, 1);
ub_TMS = zeros(numRelays, 1) + TMS_max;
lb_lambda = zeros(numRelays, 1);
ub_lambda = PS_max_i / 0.9;
lb_xi = zeros(numRelays, 1);
ub_xi = PS_max_i;
lb = [lb_TMS; lb_lambda; lb_xi];
ub = [ub_TMS; ub_lambda; ub_xi];
x0 = (lb + ub) / 2; % Initial guess in the middle of the boundary
% Since there are no linear constraints, set A, b, Aeq, beq to zero.
A = [];
bb = [];
Aeq = [];
beq = [];
%%3. Calling the fmincon solver
% We use optimoptions to select the 'interior-point' algorithm and
% set other preferences
options = optimoptions('fmincon', 'Algorithm', 'interior-point', ...
                       'Display', 'iter-detailed', ...
                       'UseParallel', false); % Utilize CPU cores if available
% Objective function and constraint definition
% Note that we use the anonymous function (@) to pass additional parameters
% (such as a, b, c, CTI_min, etc.) to the objective function and constraints.
objectiveFunction = @(x) objective_fun(x, R_pairs, numRelays, a, b, c, V_i, I_f_i, PS_max_i);
nonlinearConstraints = @(x) nonlin_con(x, R_pairs, CTI_min, t_min, numRelays, a, b, c, V_i, I_f_i, PS_max_i);
% Calling fmincon
[x_optimized, fval, exitflag, output] = fmincon(objectiveFunction, x0, A, bb, Aeq, beq, lb, ub, nonlinearConstraints, options);

total_time = 0.1168

Unrecognized function or variable 'nonlin_con'.

Error in solution>@(x)nonlin_con(x,R_pairs,CTI_min,t_min,numRelays,a,b,c,V_i,I_f_i,PS_max_i) (line 60)
nonlinearConstraints = @(x) nonlin_con(x, R_pairs, CTI_min, t_min, numRelays, a, b, c, V_i, I_f_i, PS_max_i);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in constrfunEvaluator (line 5)
[cin, ceq] = feval(Nonlcon, x, self.FunArgs.AdditionalParameters{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in OptimFunctions/constraintFirstEval (line 765)
[cin,ceq,cingrad,ceqgrad] = self.ConstraintFunAndGrad(self,self.ConFcn{3},...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fmincon (line 509)
[initVals.ncineq,initVals.nceq,initVals.gnc,initVals.gnceq,funObj] = funObj.constraintFirstEval( ...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

%% 4. Display the result
if exitflag > 0
    fprintf('\nOptimization is successfully finished!\n');
    
    % Extract the result
    optimized_TMS = x_optimized(1:numRelays);
    optimized_lambda = x_optimized(numRelays+1 : 2*numRelays);
    optimized_xi = x_optimized(2*numRelays+1 : 3*numRelays);
    
    fprintf('The value of TMS optimized: \n');
    disp(optimized_TMS);
    fprintf('The value of lambda optimized: \n');
    disp(optimized_lambda);
    fprintf('The value of xi optimized: \n');
    disp(optimized_xi);
    fprintf('Value of minimum objective function (total operation time): %.4f\n', fval);
else
    fprintf('\nOptimization unsuccessfully finished. Exit Flag: %d\n', exitflag);
    disp(output);
end
% objective_fun.m
% Objective function for fmincon.
function total_time = objective_fun(x, R_pairs, numRelays, a, b, c, V_i, I_f_i, PS_max_i)
    % Separate optimized variabel from vector x
    TMS = x(1:numRelays);
    lambda = x(numRelays+1 : 2*numRelays);
    xi = x(2*numRelays+1 : 3*numRelays);
    % Calculate PS and t for every relay
    PS = min(lambda .* V_i + xi, PS_max_i);
    I_norm = I_f_i ./ PS;    
    %  % Avoid zero dividen or negative value inside algorithm
    denominator = I_norm.^b - 1;
    t_i = TMS .* (a ./ denominator) + c;
    % Change NaN or Inf if any
    t_i(denominator <= 0) = inf;
    % Calculate total operation time from relay's pair
    total_time = 0;
    for s = 1:length(R_pairs)
        scenario_pairs = R_pairs{s};
        for k = 1:size(scenario_pairs, 1)
            i_p = scenario_pairs(k, 1); % Primary relay
            i_c = scenario_pairs(k, 2); % Backup Relay
            
            % Pastikan indeks valid
            if i_p > numRelays || i_c > numRelays
                continue; % Pass if the index is invalid
            end
            total_time = total_time + t_i(i_p) + t_i(i_c);
        end
    end
    total_time
end