summaryrefslogtreecommitdiff
path: root/code/intfmincon/ProcessSel.sce
blob: d2efefb23dcd3bb7ec15c8311c54d373ee615dea (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
//This is an MINLP example.The solution to the problem is Y1=1,Y2=0;Y3=1;C1=1;B1=1.111;B2=0;B3=1.111;BP=0;A2=0;A3=1.52;
//Netprofit = 1.923 (10^3$/hr))
//Ref:Optimization of chemical processes, second edition. By Thomas F. Edgar, David M. Himmelblau, and Leon S. Lasdon, McGraw Hill, New York, 2001, Chapter 9
//The manufacture of a chemical C in process 1 that uses raw material B.B can either be purchased or
//manufactured via two processes, 2 or 3, both of which use chemical A as a raw mate-rial.Data and 
//specifications for this example problem, involving several nonlinear input4utput relations (mass balances),
//are shown in Table.We want to deter-mine which processes to use and their production levels in order to 
//maximize profit.The processes represent design alternatives that have not yet been built. Their fixed
//costs include amortized design and construction costs over their anticipated lifetime,which are incurred
//only if the process is used
//    
//    A2------Process2-------------B2
//    A1------Process3-------------B3
//    B2+B3+BP---------------------B1
//    B1------Process1-------------C1
//    
//Problem Data
//Conversions         Process         1       C = 0.9B
//                    Process         2       B = ln(1+A)
//                    Process         3       B = 1.2ln(1+A)   (A,B,C in tons)
//                    
//Maximum Capacity    Process         1       2 ton/h of C
//                    Process         1       4 ton/h of B
//                    Process         1       5 ton/h of B
//                    
//Price       A: $1800/ton
//            B: $7000/ton
//            C: $13000/ton
//            
//Demand of C: 1 ton/h maximum
//Costs:
//---------------------------------------------------------------------------------
//                    Fixed(10^3 $/hr)        variable (10^3 $/ton of product)
//Process 1           3.5                             2
//Process 1           1                               1
//Process 1           1.5                             1.2
//---------------------------------------------------------------------------------
////======================================================================                
// Copyright (C) 2018 - IIT Bombay - FOSSEE
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// Author:Debasis Maharana
// Organization: FOSSEE, IIT Bombay
// Email: toolbox@scilab.in
//======================================================================

clc;
//Objective function to calculate profit
function costval = ProcessSel(x)
    Y1 = x(1);Y2 = x(2);Y3 = x(3);
    C1 = x(4);B1 = x(5);B2 = x(6);
    B3 = x(7);BP = x(8);A2 = x(9);A3 = x(10);

    Income = 13*C1;
    Purchase = 7 * BP;
    Chemical = 1.8*A2 + 1.8*A3;
    Invest = 3.5*Y1 + 2*C1 + Y2 + B2 + 1.5*Y3 + 1.2*B3;

    costval = -(Income - (Purchase + Chemical + Invest));
endfunction

//Nonlinear constraint function
function [C,Ceq] = Nonlincon(x)
    Y1 = x(1);Y2 = x(2);Y3 = x(3);
    C1 = x(4);B1 = x(5);B2 = x(6);
    B3 = x(7);BP = x(8);A2 = x(9);A3 = x(10);
    //No Non-linear inequality constraint
    C = [];
    Ceq(1) = C1-0.9*B1;    
    Ceq(2) = B2-log(1+A2);
    Ceq(3) = B3 - 1.2*log(1+A3);

endfunction

// Decision vector structure
//  Y1 = x(1);Y2 = x(2);Y3 = x(3);C1 = x(4);B1 = x(5);B2 = x(6);B3 = x(7);BP = x(8);A2 = x(9);A3 = x(10);

//Inequality constraints
A = [0 -4 0 0 0 1 0 0 0 0
0 0 -5 0 0 0 1 0 0 0
-2 0  0 1 0 0 0 0 0 0];

b = zeros(3,1);

//Equality Constraint
Aeq = [0 0 0 0 1 -1 -1 -1 0 0;];
beq = 0;
//Number of variables
Var = 10;
//Capacity Constraints are used for deciding upper bounds of variables
C1max = 1;
B2max = 4;
B3max = 5;

A2max = exp(B2max)-1;
A3max = exp(B3max/1.2)-1;
B1max = C1max/0.9;
BPmax = B1max;

//Bounds on variables
lb = zeros(1,Var);
// Important to give proper UB as otherwise solution might become sub optimal
//users can try ub  = [1 1 1 1 %inf %inf %inf %inf %inf %inf];

ub  = [1 1 1 C1max B1max B2max B3max BPmax A2max A3max];
//Initial guess 
x0 = lb;
//Integer ariables
intcon = [1 2 3];
mprintf('The decision variables and their bounds are")

var = ['Y1','Y2','Y3','C1','B1','B2','B3','BP','A2','A3'];
Table = [['Variables','lower Bound','Upper Bound'];[var' string(lb') string(ub')]];
disp(Table)

input("Press enter to solve the problem");
clc
//Using intfmin for solving MINLP
mprintf("Scilab is solving the problem")
[xopt,fopt,exitflag,gradient,hessian]= intfmincon(ProcessSel,x0,intcon,A,b,Aeq,beq,lb,ub,Nonlincon)
clc;
Y1 = xopt(1);Y2 = xopt(2);Y3 = xopt(3);C1 = xopt(4);B1 = xopt(5);B2 = xopt(6);
B3 = xopt(7);BP = xopt(8);A2 = xopt(9);A3 = xopt(10);

select exitflag
case 0 
    mprintf("Optimal Solution Found");
    mprintf("\n Net profit is %f $ \n Net production is : %f ton/hr",-fopt*1000,C1)
    mprintf("\n The following production process should be followed \n")
    mprintf("\n Amount of raw material A2: %f ton/hr \n Amount of raw material A3: %f ton/hr \n amount of B1 purchased : %f ton/hr",A2,A3,B1)

case 1
    mprintf("InFeasible Solution.");
case 2
    mprintf("Objective Function is Continuous Unbounded.");
case 3
    mprintf("Limit Exceeded.");
case 4
    mprintf("User Interrupt");
case 5
    mprintf("MINLP Error")
end