Need input/idea of Monte Carlo method

I am developing a system of the heat excchanger in which I am subjecting three variables to uncertainty. I labeled these in my code as p(1), p(2), and p(3). Previously, I ran fmincon to get an optimization results for this system (with those are standard variables rather than uncertainty). The feedback I received with my current method below is that I was running an optimization based on the Monte Carlo results rather than subjecting my results to a Monte Carlo analysis. I am a bit stumped here and would like some ideas or course correction to approach the method correctly:

% sample the uncertainty set and solve the optimization problem at each
% value of uncertainty
meanVal = [2.22, 3, 120]; % mean value for specific heat, flow rate, and thermal conductivity of gasoline
sigma = [0.066, 0.5, 2]; % sigma values relating to the above values

iterations = 2000; % arbitrary iteration value
OptimalL = zeros(iterations, 1); % setting up zeros to hold our data for Length [m]
OptimalR = zeros(iterations, 1); % setting up zeros to hold our data for Radius [m]

for i = 1:iterations % iterations setup
    MC = mcarlo(meanVal, sigma); % pulling monte carlo info
    LRopt = LR(MC); % getting updated values
    OptimalL(i) = LRopt(1); % instructing a new variable to take the 1st parameter (L) within the zeros we set up earlier
    OptimalR(i) = LRopt(2); % instructing a new variable to take the 2nd parameter (r) within the zeros we set up earlier
    Variable_vs(i) = LRopt(3); % instructing a new variable to take the 3rd parameter (vs) within the zeros we set up earlier
    Variable_tti(i) = LRopt(4); % instructing a new variable to take the 4th parameter (T_t(in)) within the zeros we set up earlier
    Variable_tto(i) = LRopt(5); % instructing a new variable to take the 5th parameter (T_t(out)) within the zeros we set up earlier
    Variable_tsi(i) = LRopt(6); % instructing a new variable to take the 6th parameter (T_s(in)) within the zeros we set up earlier
    Variable_tso(i) = LRopt(7); % instructing a new variable to take the 7rth parameter (T_s(out)) within the zeros we set up earlier
end

function f = obj(x) % objective function (cost of piping (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(1); % 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 = MC(3); % thermal conductivity of gasoline [W/m*K]
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 = ((x(4) - x(7)) - (x(5) - x(6))) / log((x(4) - x(7)) / (x(5) - x(6))); % deltaTLM calculation
Re_t = (rt*MC(2)*x(1))/ut; % Reynolds # of gasoline
Re_s = (rs*x(3)*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)-x(4)) - (pi*(x(1)^2)*cps)*(x(6)-x(7)) - Q; % energy balance
h(2) = x(4) - 120; % incoming gasoline will be 120 degrees C

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, [2.1534, 2.5, 118]), [2.2866, 3.5, 122]);
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;120;30;20;80]; % x = [L, r, vs, T_t(in), T_t(out), T_s(in), T_s,(out)]
options = optimoptions(@fmincon, 'Algorithm', 'SQP'); % 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

[Torsten]

And what makes you stumped ?

[Justyn Welsh]

If this were to be done correct, I would be using my results from the fmincon in a Monte Carlo analysis. Currently, I am running the Monte Carlo through fmincon, rather separately. I'm a bit confused on how to run the analysis on the results of the fmincon.

[Torsten]

If you save the fvals from the 1000 trials, you could decide after the loop which parameters constellation of the three design variables gave the best value of your objective function. Further you could plot an empirical distribution function for the fval values to see how the target variable is distributed under the assumed distributions of your uncertain input variables. Thus in contrast to @Avadhoot I think your approach is justified and correct.

Answer 1

[Avadhoot]

Hi @Justyn Welsh,

I understand that you are implementing the Monte Carlo method to check how your system performs under uncertainty in a few parameters. But the approach you are following is flawed because you are trying to introduce uncertainty before applying the optimization. But the correct approach should be the reverse of this. The feedback you received also points towards the same thing. Monte Carlo analysis is meant to subject your optimization results to uncertainty, not the other way around. The correct approach to doing this is as follows:

1. Optimization first:

• Perform optimization first without introducing uncertainty in the parameters.
• This involves setting up your "fmincon" call to use the mean values directly in your model constraints and objective function.

2. Monte Carlo analysis:

• After obtaining the optimal system you can subject it to uncertainty using the Monte Carlo method. This involves sampling the uncertain parameters (in your case, specific heat, flow rate, and thermal conductivity) from their respective distributions and evaluating how these uncertainties affect the performance or outcome of your optimized system.
• For each iteration, you sample a new set of parameters from their distributions.
• For each sampled set, evaluate how the system performs under these conditions. This does not involve re-optimizing the system; instead, you're assessing the performance of the already optimized system under uncertainty.

You need to slightly change the code to implement this approach. Before going into the iterative loop, run the "LR" function without the "MC" argument so that instead of taking the uncertain values the "fmincon" function will have all the certain values. After calculating the optimal model, then go into the iterative loop and instead of re-optimizing the model, evaluate its performance in each iteration. That way it would be the proper use of the Monte Carlo method.

I hope this helps.

[Justyn Welsh]

Hi! Thank you for the response. I have a backup code where the MC variables are not integrated into the fmincon, therefore I should be able to run the optimization process and receive the output variables.

Conceptually, I get what you're saying - that is run the optimization and then subject the response to uncertainity. I am just having a bit of trouble figuring out the MATLAB aspect of this if you might have any clues and/or hints.

For starters, is there an easy way to assign my optimization results to variables? Or is the x(1), x(2)....x(n) overwritten upon completion of the fmincon process?

Again, thank you! I look forward to figuring this out.