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)).
Monte-Carlo Analysis of fmincon optimization is converging to infeasible point
4 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)