construct a two level optimisation problem has fmincon solver for lower level and genetic algorithm for upper level

4 views (last 30 days)
Nourhan Elsayed on 16 Oct 2020
Commented: Nourhan Elsayed on 18 Oct 2020
Hi all
i have a two level optimization problem that has to run an upper level optimization problem using genetic algorithm, the upper level decision variables ""xu'' are passed to lower level as fixed parameters, for which, the lower level optimization problem should find optimimum "xl" using fmincon solver and pass these optimum "xl" to the upper level to evaluate it regarding the upper level objective function.
the code is written on 7 m.files.
to run the optimization process, run m.file.7 first, and then run m.file.4. in addition excel file for database is attached
my question is that the results showed that the upper level algorithm made only one generation and stopped with this massage.
"Optimization terminated: no feasible point found.
exitflag = -2
output = struct with fields:
problemtype: 'nonlinearconstr'
rngstate: [1×1 struct]
generations: 1
unccount: 2210
message: 'Optimization terminated: no feasible point found.'
maxconstraint: 1.2283e+03"
i am asking if my code structure is not linking the upper level and lower level on the proper way?
does the lower level pass its results to the upper level?
m.file.1 :lower level optimization solver (Low_level_opt)
function FC = Low_level_opt(xu)
% Declar Global Values
global M I J P_EL P_Prop_left P_Prop_right a b d u v w P_mcr_m P_mcr_i T IP_Lower
%% fmincon
x0 = [0.9, 0.5, 0.9, 0.4, 0.6, 0.7, 0.5, 0.7, 0.4, 0.5, 0.6, 0.7, 0.2, 0.4, 0.4, 0.5, 0.6,0.7];
lb = zeros(1,18);
ub = 0.9*ones(1,18);
options = optimoptions('fmincon', 'Algorithm' , 'sqp','PlotFcn' , 'optimplotfval')
[x,fval,exitflag,output] =fmincon(@LowerFunction,x0,[],[],[],[],lb,ub,@LLConstraint, options);
xLower = x
FC = fval
m.file.2: lower level constraints (LLConstraint)
function [c,ceq] = LLConstraint(xl)
%% Declar Global Values
global J P_mcr_G eta_Elec TEL TGEP LF PTI_left PTI_right PTO_left PTO_right P_EL
%%
% Con_Relx_F is the constraint relaxation factor and it is the limit of excess unuseful power.
Con_Relx_F = 150;
TEL = P_EL + PTI_left + PTI_right
TGEP = P_mcr_G' * LF * eta_Elec + PTO_left * eta_Elec + PTO_right * eta_Elec
c(1:J) = TEL - TGEP
c(J+1: 2*J) = TGEP - TEL -Con_Relx_F
ceq=[];
end
m.file.3 is calculating the lower level objective function (LowerFunction)
function FittingValue = LowerFunction(xl)
%% Declar Global Values
global N I a b d u v w J P_Prop_left T P_mcr_ME P_mcr_G eta_GB eta_PTO eta_PTI eta_Mech
global LF Ld_left PTI_left PTI_right PTO_left PTO_right
%% Decision Variables
Ld_left = [xl(1:6)]
LF = [xl(7:12); xl(13:18)]
%% Dependent variable
% 1) SFOC
SFOC_MLJ = a' *ones(1,J) .* Ld_left.^2 + b'*ones(1,J) .* Ld_left + d'*ones(1,J) .* ones(N,J);
SFOC_MRJ = a'*ones(1,J) .* Ld_left.^2 + b'*ones(1,J) .* Ld_left + d'*ones(1,J) .* ones(N,J);
SFOC_GIJ = u'*ones(1,J).* LF.^2 + v' *ones(1,J) .* LF + w' *ones(1,J) .* ones(I,J);
% 2) PTO and PTI
PTI_left = (P_Prop_left - Ld_left * P_mcr_ME * eta_Mech)./ (eta_Mech * eta_PTI);
PTI_left(PTI_left <0)= 0;
PTO_left = ( Ld_left * P_mcr_ME - (P_Prop_left / eta_Mech)) .* (eta_PTO * eta_GB);
PTO_left (PTO_left<0) = 0 ;
PTI_right = (P_Prop_left - Ld_left * P_mcr_ME * eta_Mech)./ (eta_Mech * eta_PTI);
PTI_right(PTI_right <0)= 0;
PTO_right = ( Ld_left * P_mcr_ME - (P_Prop_left / eta_Mech)) .* (eta_PTO * eta_GB);
PTO_right (PTO_right<0) = 0 ;
%% Objective Functions
FC = ((( P_mcr_ME * (Ld_left .* SFOC_MLJ )+ P_mcr_ME * (Ld_left .* SFOC_MRJ)) * T') + ones (1,2)* ((P_mcr_G * ones(1,J)) .* (LF .* SFOC_GIJ)) * T');
%% Output of FittingValues
FittingValue = [FC]
disp('Read of Lower fitting fn successeful')
end
the upper level problem used genetic algorithm and it also has 3 m.files
m.file.4 upper level optimization solver (opt_Up_level)
function [xu,fval]= opt_Up_level(xlower)
% Declar Global Values
global IP
PopulationSize = 100;
nvars = 3
lb = [3000, 600, 600];
ub =[10000, 6000, 6000];
options = optimoptions('ga','InitialPopulationMatrix',IP,'PopulationSize',PopulationSize, 'Display','diagnose', 'MaxGenerations',20,'PlotFcn',"gaplotbestf")
[x,fval,exitflag,output] =ga(@UpperFunction,nvars,[],[],[],[],lb,ub,@ULconstraint, [], options)
end
m.file.5 upper level constraints (ULconstraint)
function [c, ceq] = ULconstraint(xu)
% Declar Global Values
global I J P_EL P_Prop_left P_mcr_ME P_mcr_G eta_Mech eta_Elec
%% this constraints set the lower limit for Engine and Generator power
Limit_P_excess = 500;
Tot_Power = P_mcr_ME * eta_Mech + ones(1,I)* P_mcr_G * eta_Elec ;
c(2:1+J) = Tot_Power - max_Load - Limit_P_excess
ceq =[];
disp('Read of Upper Constraint is successeful')
end
m.file.6 upper level objective function (ULfunction)
function TC = UpperFunction(xu)
%% Declar Global Values
global P_mcr_G_Target P_mcr_ME_Target N I a b d u v w P_mcr_ME P_mcr_G MainEngineData_FullDatabase GeneratorData
global Fuel_Price eta_G epsilon_Conv n_Conv epsilon_MG epsilon_GB epsilon_MG epsilon_GB
%% Upper Level Decision Variables
P_mcr_ME_Target = xu(1)
P_mcr_G_Target = [xu(2:3)]
%% range of suitable Engines
P_mcr_min = 1000;
P_mcr_max = 10000;
SizeOfDatabase = size(MainEngineData_FullDatabase);
p = SizeOfDatabase(1);
z = 1;
for n = 1 : p
if P_mcr_min <= MainEngineData_FullDatabase(n,1) && MainEngineData_FullDatabase(n,1) <= P_mcr_max
MainEngineData(z,:) = MainEngineData_FullDatabase(n,:);
z = z + 1;
end
end
% Number of available engine models in the database
SizeOfDatabase = size(MainEngineData);
k = SizeOfDatabase(1);
%% Choose P_mcr_G, from generators database, closest to P_mcr_G_Target
% step 1: Writes all P_GI values from the databse in an array
P_GI_Array = GeneratorData(:,1);
% step 2: make a for-loop to find the index of P_mcr_i in the array that is close to the P_mcr_i_Target value given by the algorithm
P_mcr_G_ClosestValue = zeros(I,1);
GeneratorModel_Indices = zeros(I,1);
for n = 1 : I
P_mcr_G_ClosestValue(n) = P_GI_Array(abs(P_GI_Array-P_mcr_G_Target(n))==min(abs(P_GI_Array-P_mcr_G_Target(n))))
% Tis identifies the line of the closest value which was found before
[Value, Line]=min(abs(P_GI_Array(:)-P_mcr_G_ClosestValue(n)));
% Writes the identiefied Line in the value MainEngineModel_Indices
GeneratorModel_Indices(n) = Line
end
%% Choose P_mcr_ME ,from Engine Database, closest to target value selected by the algorithm
% step 1: Write all P_mcr_ME values from the databse in an array
P_ME_Array = MainEngineData(:,1);
% step 2: make a for-loop to find the index of P_mcr_ME in the array that is close to the P_mcr_ME_Target value given by the algorithm
P_mcr_ME_ClosestValue = zeros(N,1);
MainEngineModel_Indices = zeros(N,1);
for n = 1 : N
P_mcr_ME_ClosestValue(n) = P_ME_Array(abs(P_ME_Array-P_mcr_ME_Target(n))==min(abs(P_ME_Array-P_mcr_ME_Target(n))))
% Tis identifies the line of the closest value which was found before
[Value, Line]=min(abs(P_ME_Array(:)-P_mcr_ME_ClosestValue(n)))
% Writes the identiefied Line in the value MainEngineModel_Indices
MainEngineModel_Indices(n) = Line
end
%% Define Engine Parameters
for n = 1 : N
P_mcr_ME(n) = MainEngineData(MainEngineModel_Indices(n),1)
a(n) = MainEngineData(MainEngineModel_Indices(n),3);
b(n) = MainEngineData(MainEngineModel_Indices(n),4);
d(n) = MainEngineData(MainEngineModel_Indices(n),5);
end
% Define Generators Parameter
for n = 1 : I
P_mcr_G(n) = GeneratorData(GeneratorModel_Indices(n),1)
u(n) = GeneratorData(GeneratorModel_Indices(n),3);
v(n) = GeneratorData(GeneratorModel_Indices(n),4);
w(n) = GeneratorData(GeneratorModel_Indices(n),5)
end
%% Objective Function
IC = (250* P_mcr_ME + 2000000) + ones(1,I) * (800 * ((P_mcr_G / 1000).^-6) .* (P_mcr_G / eta_G)) + epsilon_Conv * n_Conv + epsilon_MG * 2 + epsilon_GB * 2
FC = @Low_level_opt;
TFC = FC(xu)
OC = TFC * Fuel_Price
TC = OC + IC
disp('Read of Upper fitting fn successeful')
end
m.file.7 containes input data and it is linked to excel file (InputData)
function [P_mcr_G_Target, P_mcr_ME_Target, N, I, IP_Lower, J, P_EL, P_Prop_left, P_Prop_right, T, eta_SWBD, P_mcr_ME, P_mcr_G, eta_H, eta_B, eta_S, eta_GB, Fuel_Price, eta_G, eta_PTO, eta_PTI, eta_Conv, Vs_J, P_ELSafelimit, epsilon_Conv, n_Conv, epsilon_MG, epsilon_GB, Ld_left, LF, eta_Mech, eta_Elec, PTI_left, PTI_right, PTO_left, PTO_right, SFOC_MLJ, SFOC_MRJ, SFOC_GIJ] = InputData(N, I, ShipData, MainEngineData_FullDatabase, GeneratorData, OtherInputData)
%% Declar Global Values
global P_mcr_G_Target P_mcr_ME_Target N I IP_Lower
global J P_EL P_Prop_left P_Prop_right T IP P_Prop_Max P_EL_Max eta_SWBD P_mcr_ME P_mcr_G
global eta_H eta_B eta_S eta_GB ShipData MainEngineData_FullDatabase GeneratorData OtherInputData
global Fuel_Price eta_G eta_PTO eta_PTI eta_Conv Vs_J P_ELSafelimit epsilon_Conv n_Conv epsilon_MG epsilon_GB
global epsilon_MG epsilon_GB Ld_left LF
global eta_Mech eta_Elec PTI_left PTI_right PTO_left PTO_right SFOC_MLJ SFOC_MRJ SFOC_GIJ
%% Data Entery
N = 1;
I = 2;
% ship input data (Profil)
T = ShipData(5,:);
P_Prop_left = ShipData(6,:);
P_Prop_right = ShipData(7,:);
P_EL = ShipData(8,:);
P_Prop_Max = max(P_Prop_left);
P_EL_Max = max(P_EL);
SizeOfOperationalStages = size(P_Prop_left);
J = SizeOfOperationalStages(2);
eta_G = OtherInputData(6);
eta_H = OtherInputData(7);
eta_B = OtherInputData(8);
eta_S = OtherInputData(9);
eta_GB = OtherInputData(10);
eta_SWBD = OtherInputData(11);
eta_PTO = OtherInputData(12);
eta_PTI = OtherInputData(13);
eta_Conv = OtherInputData(14);
Fuel_Price = OtherInputData(15);
Vs_J = ShipData(3,:);
P_ELSafelimit = OtherInputData(1);
epsilon_Conv= OtherInputData(2);
n_Conv = OtherInputData(3);
epsilon_MG = OtherInputData(4);
epsilon_GB = OtherInputData(5);
%% Total efficiency of propulsion train
eta_Mech = eta_B * eta_H * eta_S * eta_GB;
eta_Elec = eta_SWBD * eta_Conv;
%% IP
StD1 = 1000;
StD2 = 500;
P_mcr_ME_Target_initial = (max(P_Prop_left)*1.2) / N
P_mcr_G_Target_initial = (P_EL_Max * 1.1)/I
Population_Size = 50
Pop_Matrix = [P_mcr_ME_Target_initial, P_mcr_G_Target_initial]
IP_sec = repmat(Pop_Matrix, [Population_Size,1]);
if N==1 && I== 2
IP_sec(2:end,1) = IP_sec(2:end,1) + StD1*randn(Population_Size-1,1)
IP_sec(2:end,2) = IP_sec(2:end,2) + StD2*randn(Population_Size-1,1)
IP = [IP_sec, IP_sec(:,2)];
IP(IP<0) = (IP(IP<0) * -1)
end
if N==2 && I==2
IP = [IP_sec(:,1), IP_sec(:,1), IP_sec(:,2), IP_sec(:,2)];
IP(IP<0) = (IP(IP<0) * -1)
end
%% define variables structure
SFOC_MLJ = zeros(N,J);
SFOC_MRJ = zeros(N,J);
SFOC_GIJ = zeros(I,J);
PTI_left = ones(1,J);
PTI_right = ones(1,J);
PTO_left= ones(1,J);
PTO_right= ones(1,J);
P_mcr_ME = ones(N,1);
P_mcr_G = ones(I,1);
P_mcr_ME_Target = ones(N, 1);
P_mcr_G_Target = ones(I, 1);
Ld_left = ones(N,J);
LF = ones(I,J);
end

Matt J on 18 Oct 2020
Using global variables to pass around constant data is quite inadvisable. You should use anonymous or nested functions instead.
Nourhan Elsayed on 18 Oct 2020
Thank you very much for this link

Matt J on 18 Oct 2020
Edited: Matt J on 18 Oct 2020
You are passing a variable called xl to LLConstraint(), but xl is not used anywhere within that function. The function is therefore making a judgement about whether xl is feasible or not completely independently of xl.

Nourhan Elsayed on 18 Oct 2020
so, what you suggest??
Matt J on 18 Oct 2020
Get rid of all your xl-dependent global variables. Explicitly compute any xl-dependent quantities within each separate function that needs them, and don't try to share the results between functions by making them global.
Nourhan Elsayed on 18 Oct 2020
OK. i will make your proposal and back to tell you the result
thank you Matt J