Monte-Carlo Analysis of fmincon optimization is converging to infeasible point
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Justyn Welsh
el 29 de Feb. de 2024
Respondida: Torsten
el 29 de Feb. de 2024
I put together an fmincon optimization code that worked fine and properly convergenced. I wanted to add uncertainty into this, however; it no longer converges to a feasible point. I'm confused as when I have ran the numbers involved with the uncertainty in the original optimization program, it would converge each time. Looking for potential causes of this issue:
% sample the uncertainty set and solve the optimization problem at each
% value of uncertainty
meanVal = [125; 2.22]; % mean value for temperature of incoming gasoline & specific heat of gasoline
sigma = [1; 0.03]; % sigma values relating to the above values
iterations = 2; % arbitrary iteration value
OptimalLR = zeros(iterations, 1); % setting up zeros to hold our data
for i = 1:iterations % iterations setup
MC = mcarlo(meanVal, sigma); % pulling monte carlo info
LRopt = LR(MC); % getting updated values
OptimalLR(i) = LRopt(1); % instructing a new variable to take the 6th parameter (V) within the zeros we set up earlier
end
function f = obj(x) % objective function (cost pipiing (radius and length) to be minimized)
f = ((0.10*x(2))+1)*(15*x(1));
end
function [g,h] = model(x,MC)
% Implement the constraints that define the model
% Use h for equality constraints and g for inequality constraints
cpt = MC(2); % specific heat of gasoline [kJ/kgK] - MC(2) is specific heat of gasoline uncertainity
cps = 4.184; % specific heat of water [kJ/kgK]
rt = 730; % density of gasoline [kg/m^3]
rs = 997; % density of water [kg/m^3]
ut = 20; % dynamic viscosity of gasoline [mPa*s] - this value is at roughly 120 C for gasoline (uncertainity spans this range)
us = 1; % dynamic viscosity of water [mPa*s] - this value is at the temperature of my initial guess for water
kt = 120; % thermal conductivity of gasoline [W/m*K] - maybe uncertainty?
ks = 0.598; % thermal conductivity of water [W/m*K]
h_f = 0.001; % thermal resistance of the fouling [K/W] - provided constant based off average fouling value for stainless steel heat exchangers
% all constant values
t_LM = ((MC(1) - x(7)) - (x(5) - x(6))) / log((MC(1) - x(7)) / (x(5) - x(6))); % deltaTLM calcuation - MC(1) is temperature incoming of gasoline uncertainty
Re_t = (rt*x(3)*x(1))/ut; % Reynolds # of gasoline
Re_s = (rs*x(4)*x(1))/us; % Reynolds # of water
Pr_t = (ut*cpt)/kt; % Prandtl # of gasoline
Pr_s = (us*cps)/ks; % Prandtl # of water
Nu_t = 0.023*(Re_t^0.8)*(Pr_t^0.4); % Nusselt # of gasoline
Nu_s = 0.023*(Re_s^0.8)*(Pr_s^0.4); % Nusselt # of water
h_t = (Nu_t*kt)/x(1); % thermal coefficient of gasoline
h_s = (Nu_s*ks)/x(1); % thermal coefficient of water
R_t = 1/(h_t*pi*(x(2)^2)); % thermal resistance of the tube
R_s = 1/(h_s*pi*(x(2)^2)); % thermal resistance of the shell
R_f = 1/(h_f*pi*(x(2)^2)); % thermal resistance of the fouling
R_total = 1/((1/R_t) + (1/R_s) + (1/R_f)); % total thermal resistance
U = 1/R_total; % overall heat transfer coefficient
Q = U*pi*(x(2)^2)*t_LM; % heat value
p = [20, 22, 150]; % parameters for nonlinear constraints
% all relevant equations to the system
h = zeros(2,1); % equality constraints
h(1) = (pi*(x(1)^2)*cpt)*(x(5)-MC(1)) - (pi*(x(1)^2)*cps)*(x(6)-x(7)) - Q; % energy balance
g = zeros(3,1); % inequality constraints - p is referenced in fmincon setup
g(1) = p(1) - x(5); % final temperature of gasoline cannot exceed 22
g(2) = x(5) - p(2); % final temperature of gasoline must exceed 20
g(3) = p(3) - ((0.10*x(2))+1)*(15*x(1)); % total cost of piping cannot exceed $150 based on piping length being $15/m and radius adding 10% to the cost to scale
end
function MC = mcarlo(meanVal, sigma) % setting up monte carlo method
MC = meanVal + sigma.*randn(size(meanVal)); % inspired from p=sigma.*randn+kmean from lecture
MC = min(max(MC, [120, 2.1534]), [130, 2.2866]);
end
function LRopt = LR(MC) % this becomes our new optimization problem now that we are subjecting our variables to uncertainty
% solve the optimization problem from here
f = @(x)obj(x); % declare objective function
lb = [0, 0, -Inf, -Inf, -Inf, 0, 0]; % lower bounds on x
ub = []; % upper bounds on x
A = []; % linear inequality constraints - NONE
b = []; % NONE
Aeq = []; % linear equality constraints - NONE
beq = []; % NONE
nonlcon = @(x)model(x,MC); % model for nonlinear constraints
% initial guess and algorithm
x0 = [10;2;3;3;30;20;80]; % x = (L, r, vt, vs, T_t(out), T_s(in), T_s,(out))
options = optimoptions(@fmincon, 'Algorithm', 'SQP', 'Display','iter'); % use SQP algorithm
% call the SQP solver to find solution to the problem
LRopt = fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
end
0 comentarios
Respuesta aceptada
Torsten
el 29 de Feb. de 2024
I guess log((MC(1) - x(7)) becomes complex-valued because MC(1)-x(7) becomes negative. Complex-valued constraints have to be avoided (e.g. by putting a constraint on x(7) that it has to be smaller than MC(1)).
0 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!