__Newton’s Method using Matlab, Chemical Reaction Engineering__

This sample matlab assignment solution showcases the derivation of newton’s method using matlab. It has been used to solve problems in chemical reaction engineering subjects. In this example, two stirred tank reactors are connected in series. One compound is being converted into another in two reactors. The unconverted compound and the product are being separated in a separation unit and the unconverted compound is being recycled back. The expert states that the total cost of the system is a function of conversion factors in two reactors. The tutor has demonstrated to find out optimal conversion factors in two reactors for which the total cost of the system is minimum. The cost function in terms of two conversion factors is given as follows,

The expert has derived a version of Newton’s method for solving this problem and has compared it with the solution obtained by fminsearch. He has also calculated the Hessian.

** SOLUTION** : –

clc

clear

format long

% Function Definition (Enter your Function here):

syms X Y;

f = (X/(Y*(1+X)^2))^0.6+((1-X/Y)/(1-Y)^2)^0.6+6*(1/Y)^0.6;

% Initial Guess (Choose Initial Guesses):

x(1) = 0.2;

y(1) = 0.4;

e = 10^(-8); % Convergence Criteria

i = 1; % Iteration Counter

% Gradient and Hessian Computation:

df_dx = diff(f, X);

df_dy = diff(f, Y);

J = [subs(df_dx,[X,Y], [x(1),y(1)]) subs(df_dy, [X,Y], [x(1),y(1)])]; % Gradient

ddf_ddx = diff(df_dx,X);

ddf_ddy = diff(df_dy,Y);

ddf_dxdy = diff(df_dx,Y);

ddf_ddx_1 = subs(ddf_ddx, [X,Y], [x(1),y(1)]);

ddf_ddy_1 = subs(ddf_ddy, [X,Y], [x(1),y(1)]);

ddf_dxdy_1 = subs(ddf_dxdy, [X,Y], [x(1),y(1)]);

H = [ddf_ddx_1, ddf_dxdy_1; ddf_dxdy_1, ddf_ddy_1]; % Hessian

S = inv(H); % Search Direction

% Optimization Condition:

while norm(J) > e

I = [x(i),y(i)]’;

x(i+1) = I(1)-S(1,:)*J’;

y(i+1) = I(2)-S(2,:)*J’;

i = i+1;

J = [subs(df_dx,[X,Y], [x(i),y(i)]) subs(df_dy, [X,Y], [x(i),y(i)])]; % Updated Jacobian

ddf_ddx_1 = subs(ddf_ddx, [X,Y], [x(i),y(i)]);

ddf_ddy_1 = subs(ddf_ddy, [X,Y], [x(i),y(i)]);

ddf_dxdy_1 = subs(ddf_dxdy, [X,Y], [x(i),y(i)]);

H = [ddf_ddx_1, ddf_dxdy_1; ddf_dxdy_1, ddf_ddy_1]; % Updated Hessian

S = inv(H); % New Search Direction

end

% Result Table:`

Iter = 1:i;

X_coordinate = x’;

Y_coordinate = y’;

Iterations = Iter’;

T = table(Iterations,X_coordinate,Y_coordinate);

% Plots:

fcontour(f, ‘Fill’, ‘On’);

hold on;

plot(x,y,’*-r’);

grid on;

% Output:

fprintf(‘Initial Objective Function Value: %d\n\n’,subs(f,[X,Y], [x(1),y(1)]));

if (norm(J) < e)

fprintf(‘Minimum succesfully obtained…\n\n’);

end

fprintf(‘Number of Iterations for Convergence: %d\n\n’, i);

fprintf(‘Point of Minima: [%d,%d]\n\n’, x(i), y(i));

fprintf(‘Objective Function Minimum Value after Optimization: %f\n\n’, subs(f,[X,Y], [x(i),y(i)]));

disp(T)

% Cost = @(X,Y) (X./(Y.*(1+X)^2)).^0.6+((1-X./Y)/(1-Y).^2)^0.6+6*(1/Y).^0.6;

v=[0.2,0.4];

w=fminsearch(@two_var,v);

function b = two_var(v)

x = v(1);

y = v(2);

b = (x./(y.*(1+x)^2)).^0.6+((1-x./y)/(1-y).^2)^0.6+6*(1/y).^0.6;