## 65  Methanol to Hydrocarbons

Benchmarking Optimization Software with COPS Elizabeth D. Dolan and Jorge J. More ARGONNE NATIONAL LABORATORY

### 65.1  Problem Formulation

Find theta over t in [0; 1.122] to minimize

J =
 3 ∑ j=1

 21 ∑ i=1
(yj,i − yj,i,meas)2

subject to:

 dy1 dt
= −(2*theta2
 theta1*y2 (theta2+theta5)*y1+y2
+theta3+theta4)*y1
 dy2 dt
=
 theta1*y1*(theta2*y1−y2) (theta2+theta5)*y1+y2
+theta3*y1
 dy3 dt
=
 theta1*y1*(y2+theta5*y1) (theta2+theta5)*y1+y2
+theta4*y1
 theta >= 0

Where the data is given in the code.

Reference: [14]

### 65.2  Problem setup

```toms t theta1 theta2 theta3 theta4 theta5

% Various constants and expressions
y1meas = [0.7085;0.5971;0.5537;0.3684;0.1712;...
0.1198;0.0747;0.0529;0.0415;0.0261;0.0208;...
0.0085;0.0053;0.0019;0.0018];
y2meas = [0.1621;0.1855;0.1989;0.2845;0.3491;...
0.3098;0.3576;0.3347;0.3388;0.3557;0.3483;...
0.3836;0.3611;0.3609;0.3485];
y3meas = [0.0811;0.0965;0.1198;0.1535;0.2097;...
0.2628;0.2467;0.2884;0.2757;0.3167;0.2954;...
0.295;0.2937;0.2831;0.2846];
tmeas = [0.05;0.065;0.08;0.123;0.233;0.273;...
0.354;0.397;0.418;0.502;0.553;...
0.681;0.75;0.916;0.937];
```

### 65.3  Solve the problem, using a successively larger number collocation points

```for n=[20 80]
```
```    p = tomPhase('p', t, 0, 1.122, n);
setPhase(p);

tomStates y1 y2 y3

% Initial guess
if n == 20
x0 = {theta1 == 1; theta2 == 1
theta3 == 1; theta4 == 1
theta5 == 1
icollocate({
y1 == 1-(1-0.0006)*t/1.122
y2 == 0.3698*t/1.122
y3 == 0.2899*t/1.122})};
else
x0 = {theta1 == theta1opt; theta2 == theta2opt
theta3 == theta3opt; theta4 == theta4opt
theta5 == theta5opt
icollocate({
y1 == y1opt
y2 == y2opt
y3 == y3opt})};
end

% Box constraints
cbox = {sqrt(eps) <= theta1; sqrt(eps) <= theta2
sqrt(eps) <= theta3; sqrt(eps) <= theta4
sqrt(eps) <= theta5};

y1err = sum((atPoints(tmeas,y1) - y1meas).^2);
y2err = sum((atPoints(tmeas,y2) - y2meas).^2);
y3err = sum((atPoints(tmeas,y3) - y3meas).^2);

% Start and end points cannot be interpolated
y1end = (1-initial(y1)).^2 + (0.0006-final(y1))^2;
y2end = (0-initial(y2)).^2 + (0.3698-final(y2))^2;
y3end = (0-initial(y3)).^2 + (0.2899-final(y3))^2;

% ODEs and path constraints
ceq = collocate({
dot(y1) == -(2*theta2-(theta1*y2)./((theta2+theta5)*y1+y2)+theta3+theta4).*y1
dot(y2) == (theta1*y1.*(theta2*y1-y2))./((theta2+theta5)*y1+y2)+theta3*y1
dot(y3) == (theta1*y1.*(y2+theta5*y1))./((theta2+theta5)*y1+y2)+theta4*y1});

% Objective
objective = y1err+y2err+y3err+y1end+y2end+y3end;
```

### 65.4  Solve the problem

```    options = struct;
options.name = 'Methanol to Hydrocarbons';
solution = ezsolve(objective, {cbox, ceq}, x0, options);

% Optimal x, theta for starting point
y1opt = subs(y1, solution);
y2opt = subs(y2, solution);
y3opt = subs(y3, solution);
theta1opt = subs(theta1, solution);
theta2opt = subs(theta2, solution);
theta3opt = subs(theta3, solution);
theta4opt = subs(theta4, solution);
theta5opt = subs(theta5, solution);
```
```Problem type appears to be: qpcon
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - Tomlab Optimization Inc. Development license  999001. Valid to 2011-02-05
=====================================================================================
Problem: ---  1: Methanol to Hydrocarbons       f_k       0.008301664004164877
sum(|constr|)      0.000000001050742318
f(x_k) + sum(|constr|)      0.008301665054907195
f(x_0)     -0.959232294294469990

Solver: snopt.  EXIT=0.  INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv    1 ConstrEv   42 ConJacEv   42 Iter   41 MinorIter   88
CPU time: 0.140625 sec. Elapsed time: 0.157000 sec.
```
```Problem type appears to be: qpcon
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - Tomlab Optimization Inc. Development license  999001. Valid to 2011-02-05
=====================================================================================
Problem: ---  1: Methanol to Hydrocarbons       f_k       0.008301664004168430
sum(|constr|)      0.000000988210502411
f(x_k) + sum(|constr|)      0.008302652214670841
f(x_0)     -5.007954925995837100

Solver: snopt.  EXIT=0.  INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv    1 ConstrEv    1 ConJacEv    1 MinorIter  159
CPU time: 0.125000 sec. Elapsed time: 0.125000 sec.
```
```end

t  = subs(collocate(t),solution);
y1 = collocate(y1opt);
y2 = collocate(y2opt);
y3 = collocate(y3opt);
t1 = subs(theta1,solution);
t2 = subs(theta2,solution);
t3 = subs(theta3,solution);
t4 = subs(theta4,solution);
t5 = subs(theta5,solution);
```

### 65.5  Plot result

```figure(1);
tm  = [0;tmeas;1.122];
y1m = [1;y1meas;0.0006];
y2m = [0;y2meas;0.3698];
y3m = [0;y3meas;0.2899];
plot(t,y1,'*-',t,y2,'*-',t,y3,'*-',tm,y1m,'ko',tm,y2m,'ko',tm,y3m,'ko');
legend('y1','y2','y3','ymeas');
title(sprintf('Methanol to Hyd, theta = [%g %g %g %g %g]',t1,t2,t3,t4,t5));
```