summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--code/Symphonymat/Factory planing.sce143
-rw-r--r--code/Symphonymat/Mining.sce154
-rw-r--r--code/Symphonymat/factory.xlsbin0 -> 29184 bytes
-rw-r--r--code/fminbnd/Chichinadze.sce64
-rw-r--r--code/fminbnd/Cola.sce89
-rw-r--r--code/fmincon/Tankprob.sci95
-rw-r--r--code/fmincon/spring.sci120
-rw-r--r--code/fminsearch/Brownsfunc.sci73
-rw-r--r--code/fminsearch/FletcherPowell.sce65
-rw-r--r--code/fminunc/Freudenstein and Roth.sci75
-rw-r--r--code/fsolve/Redlich-Kwong.sce69
-rw-r--r--code/fsolve/SLEs.sce50
-rw-r--r--code/intfmincon/ProcessSel.sce148
-rw-r--r--code/intqpipopt/SpringCart.sce68
-rw-r--r--code/linprog/Blending_problem.sce135
-rw-r--r--code/linprog/Employee Scheduling.sce83
-rw-r--r--code/linprog/FuelOil.sce86
-rw-r--r--code/linprog/Reddy_Mikks.sce84
-rw-r--r--code/linprog/TOOLCO.sce82
-rw-r--r--code/linprog/TOYCO Dual.sce88
-rw-r--r--code/linprog/TOYCO model.sce67
21 files changed, 1838 insertions, 0 deletions
diff --git a/code/Symphonymat/Factory planing.sce b/code/Symphonymat/Factory planing.sce
new file mode 100644
index 0000000..e140f42
--- /dev/null
+++ b/code/Symphonymat/Factory planing.sce
@@ -0,0 +1,143 @@
+//A practical problem of factory planning which determines the optimum product
+// mix subject to production capacity and marketing limitation.
+//This example shows how to use spreadsheet data directly in scilab.
+//Ref: H. Paul Williams ,Model Building in Mathematical Programming ,
+//A John Wiley & Sons, Ltd., Publication, Fifth Ed, Chapter 12.
+//
+//Example: An engineering factory makes seven products (PROD 1 to PROD 7) on the
+//following machines: four grinders, two vertical drills, three horizontal drills, one
+//borer and one planer. Each product yields a certain contribution to profit (defined
+//as £/unit selling price minus cost of raw materials). These quantities (in £/unit)
+//together with the unit production times (hours) required on each process are given
+//below. A dash indicates that a product does not require a process.
+//---------------------------------------------------------------------------------------
+// Prod 1 Prod 2 Prod 3 Prod 4 Prod 5 Prod 6 Prod 7
+//
+//Contribution to profit 10 6 8 4 11 9 3
+//Grinding 0.5 0.7 - - 0.3 0.2 0.5
+//Vertical drilling 0.1 0.2 - 0.3 - 0.6 -
+//Horizontal drilling 0.2 - 0.8 - - - 0.6
+//Boring 0.05 0.03 - 0.07 0.1 - 0.08
+//Planing - - 0.01 - 0.05 - 0.05
+//--------------------------------------------------------------------------------------
+//In the present month (January) and the five subsequent months, certain
+//machines will be down for maintenance. These machines will be as follows:
+//--------------------------------------------
+//January 1 Grinder
+//February 2 Horizontal drills
+//March 1 Borer
+//April 1 Vertical drill
+//May 1 Grinder and 1 Vertical drill
+//June 1 Planer and 1 Horizontal drill
+//--------------------------------------------
+//
+//There are marketing limitations on each product in each month. These are
+//given in the following table:
+//-------------------------------------------
+// 1 2 3 4 5 6 7
+//-------------------------------------------
+//January 500 1000 300 300 800 200 100
+//February 600 500 200 0 400 300 150
+//March 300 600 0 0 500 400 100
+//April 200 300 400 500 200 0 100
+//May 0 100 500 100 1000 300 0
+//June 500 500 100 300 1100 500 60
+//-------------------------------------------
+//
+//The factory works at six days a week with two shifts of 8 h each day.
+//When and what should the factory make in order to maximise the total profit?
+//This problem considers single period only with no storage
+
+// 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
+
+filepath = 'C:\Users\Iball\Desktop\scilab problems\Date 21-05-2018-Debasis\factory.xls';
+[fd,SST,Sheetnames,Sheetpos] = xls_open(filepath);
+S = readxls(filepath);
+[Contributation,TextInd] = xls_read(fd,Sheetpos(Sheetnames == 'Contributation')); //Sheetpos gives the position of all the sheets
+[Available,TextInd] = xls_read(fd,Sheetpos(Sheetnames == 'Available'));
+[Failure,TextInd] = xls_read(fd,Sheetpos(Sheetnames == 'Machine Failure'));
+[Limitation,TextInd] = xls_read(fd,Sheetpos(Sheetnames == 'Limitation'));
+mclose(fd) //close the file
+
+Sheetdata = readxls(filepath);
+mprintf('Problem Data received\n')
+for i = 1:4
+ disp(Sheetnames(i))
+ disp(Sheetdata(i))
+end
+input('press enter to continue')
+clc
+mprintf('Scilab is solving your problem')
+Shift_T = 8;N_shift = 2;N_Workingdys = 24;//Data is fixed according to example. However it can be made as user defined parameter
+
+Wrk_hrs = Shift_T*N_shift*N_Workingdys;
+
+[Nmonth,Nprod] = size(Limitation);
+
+Nprod = Nprod - 1;
+Nmonth = Nmonth - 1;
+Ntype = size(Available,'r');
+
+Available = Available(:,2);
+
+Contributation = Contributation(2:Ntype+2,2:Nprod+1);// Extracting the data matrix only
+Limitation = Limitation(2:Nmonth+1,2:Nprod+1);
+Failure = Failure(2:Nmonth+1,2:Ntype+1);
+Failure(isnan(Failure))=0;
+
+ub = [];
+for i = 1:Ntype
+ Amonth(i,:) = Contributation(i+1,:);
+ breq(i) = Available(i)*Wrk_hrs;
+end
+b = [];A = zeros(Nmonth*Ntype,Nmonth*Nprod);C = [];lb = zeros(1,Nmonth*Nprod);
+
+// Instead of seperate loops, we can use a single loop for determining all the possible variables
+for i = 1:Nmonth
+ A((i-1)*Ntype+1:i*Ntype,(i-1)*Nprod+1:i*Nprod) = Amonth;
+ B_month = breq - (Failure(i,:))'*Wrk_hrs;
+ b = [b;B_month];
+ C = [C;-Contributation(1,:)'];
+ ub = [ub Limitation(i,:)];
+end
+
+intcon = 1:Nmonth*Nprod; // all production units are integers
+options = list("time_limit", 2500);
+[xopt,fopt,status,output] = symphonymat(C,intcon,A,b,[],[],lb,ub,options);
+clc
+select status
+case 227 then
+ mprintf('Optimal Solution Found')
+case 228 then
+ mprintf('Maximum CPU Time exceeded')
+case 229 then
+ mprintf('Maximum Number of Node Limit Exceeded')
+else
+ mprintf('Maximum Number of Iterations Limit Exceeded')
+end
+input('Press enter to view results')
+// Solution representation
+A2 = ['January','Februry','March','April','May','June']';// users can modify it to accept from the excel sheet
+A1 = [" ", 'Prod1','Prod2','Prod3','Prod4','Prod5','Prod6','Prod7'];
+
+for i = 1:Nmonth
+ solution(i,:) = string(xopt((i-1)*Nprod+1:i*Nprod)');
+end
+
+mprintf('Production schedule for all the months')
+
+table = [A1;[A2 solution]];
+disp(table)
+mprintf('The profit is %d',-fopt)
+
+
diff --git a/code/Symphonymat/Mining.sce b/code/Symphonymat/Mining.sce
new file mode 100644
index 0000000..f92b93a
--- /dev/null
+++ b/code/Symphonymat/Mining.sce
@@ -0,0 +1,154 @@
+// This is an example of mixed integer linear programming problem. This problem dimension increases significantly with increase in years or mines.
+//Example: A mining company is going to continue operating in a certain area for the next five years. There are four mines in this area, but it can operate at most three in any one year. Although a mine may not operate in a certain year, it is still necessary to keep it ‘open’, in the sense that royalties are payable, if it be operated in a future year. Clearly, if a mine is not going to be worked again, it can be permanently closed down and no moreroyalties need be paid. The yearly royalties payable on each mine kept ‘open’ are as follows:
+// -----------------
+// Mine 1 £5 million
+// Mine 2 £4 million
+// Mine 3 £4 million
+// Mine 4 £5 million
+// ------------------
+//There is an upper limit to the amount of ore, which can be extracted from each mine in a year. These upper limits are as follows:
+// ---------------------
+// Mine 1 2 × 10^6 tons
+// Mine 2 2.5 × 10^6 tons
+// Mine 3 1.3 × 10^6 tons
+// Mine 4 3 × 10^6 tons
+// ---------------------
+//The ore from the different mines is of varying quality. This quality is measured on a scale so that blending ores together results in a linear combination of the quality measurements, for example, if equal quantities of two ores were combined, the resultant ore would have a quality measurement half way between that of the ingredient ores. Measured in these units the qualities of the ores from the mines are given as follows:
+// ---------------
+// Mine 1 1.0
+// Mine 2 0.7
+// Mine 3 1.5
+// Mine 4 0.5
+// ---------------
+//In each year, it is necessary to combine the total outputs from each mine to produce a blended ore of exactly some stipulated quality. For each year, these qualities are as follows:
+// -------------------
+// Year 1 0.9
+// Year 2 0.8
+// Year 3 1.2
+// Year 4 0.6
+// Year 5 1.0
+// -------------------
+//The final blended ore sells for £10 ton each year. Which mines should be operated each year and how much should they produce?
+
+// 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
+
+Nmines = 4;
+//Royltes = [5e6 4e6 4e6 5e6];
+Roylte = [5e6 4e6 4e6 5e6];
+
+ubMines = [2e6 2.5e6 1.3e6 3e6];
+Quality = [1 0.7 1.5 0.5];
+Quality_blend = [0.9 0.8 1.2 0.6 1];
+Nyears = length(Quality_blend);
+Royltes = Roylte;
+for i = 1:Nyears-1
+ Roylte = Roylte*0.9;
+ Royltes = [Royltes Roylte];
+end
+
+
+
+mprintf('Received Data')
+mprintf('\nNumber of years to operate %d',Nyears)
+mprintf('\nNumber of mines %d',Nmines);
+disp(Royltes,'Royalities to be paid for each mine',);
+disp(ubMines,'Upper limit of ore available from each mine');
+disp(Quality,'Quality of ore from each mine');
+input('Press enter to proceed')
+
+//i - mine , t - year
+A = zeros(1,3*Nyears*Nmines+Nyears)
+Aeq_quantity = zeros(Nyears,Nyears*(3*Nmines+1));
+Aeq_quality = zeros(Nyears,Nyears*(3*Nmines+1));
+A_operation = zeros(Nyears,Nyears*(3*Nmines+1));
+Aoutput =[];Aopen=[];
+boutput = zeros(Nyears*Nmines,1);
+bopen = zeros(Nyears*Nmines,1);
+beq_quantity = zeros(Nyears,1);
+beq_quality = zeros(Nyears,1);
+
+for i = 1:Nyears
+ Aeq_quantity(i,(i-1)*Nmines+1:i*Nmines) = 1;
+ Aeq_quantity(i,Nyears*Nmines+i) = -1;
+ beq_quantity(i) = 0;
+
+ Aeq_quality(i,(i-1)*Nmines+1:i*Nmines) = Quality ;
+ Aeq_quality(i,Nmines*Nyears+i) = -Quality_blend(i);
+ beq_quality(i) = 0
+
+ A_operation(i,Nyears*(Nmines+1)+((i-1)*Nmines+1:i*Nmines)) = 1;
+ B_operation(i) = 3;
+ A_output = zeros(Nmines,Nyears*(3*Nmines+1));
+ for j = 1:Nmines
+ A_output(j,(i-1)*Nmines+j) = 1;
+ A_output(j,Nyears*(Nmines+1)+(i-1)*Nmines+j) = - ubMines(j);
+ end
+ Aoutput = [Aoutput;A_output];
+ A_open = zeros(Nmines,Nyears*(3*Nmines+1));
+ for j = 1:Nmines
+ A_open(j,(Nmines+1)*Nyears+(i-1)*Nmines+j) = 1;
+ A_open(j,(2*Nmines+1)*Nyears+(i-1)*Nmines+j) = -1;
+ end
+ Aopen = [Aopen;A_open];
+end
+Aclose = [];bclose = zeros((Nyears-1)*Nmines,1);
+for i = 1:Nyears-1
+ A_close = zeros(Nmines,Nyears*(3*Nmines+1));
+ for j = 1:Nmines
+ A_close(j,Nyears*(2*Nmines+1)+(i-1)*Nmines+j) = -1;
+ A_close(j,Nyears*(2*Nmines+1)+i*Nmines+j) = 1;
+ end
+ Aclose = [Aclose;A_close];
+end
+
+A = [Aoutput;A_operation;Aopen;Aclose]
+Aeq = [Aeq_quantity;Aeq_quality];
+b = [boutput;B_operation;bopen;bclose]
+beq = [beq_quantity;beq_quality];
+
+lb = zeros(1,Nyears*(3*Nmines+1));
+ub = [repmat(ubMines,1,Nyears) sum(ubMines)*ones(1,Nyears) ones(1,2*Nmines*Nyears)];
+//C = [zeros(1,Nmines*Nyears) [10 9 8 7 6].*ones(1,Nyears) zeros(1,Nyears*Nmines) -repmat(Royltes,1,Nyears)]';
+C = [zeros(1,Nmines*Nyears) [10 10*0.9 10*0.9^2 10*0.9^3 10*0.9^4 ].*ones(1,Nyears) zeros(1,Nyears*Nmines) -Royltes]';
+
+// All variables are not integers
+intcon = [(Nmines+1)*Nyears+1:Nyears*(3*Nmines+1)];
+options = list("time_limit", 2500);
+[xopt,fopt,status,output] = symphonymat(-C,intcon,A,b,Aeq,beq,lb,ub,options);
+clc
+select status
+case 227 then
+ mprintf('Optimal Solution Found')
+case 228 then
+ mprintf('Maximum CPU Time exceeded')
+case 229 then
+ mprintf('Maximum Number of Node Limit Exceeded')
+else
+ mprintf('Maximum Number of Iterations Limit Exceeded')
+end
+input('Press enter to view results')
+
+for i = 1:Nyears
+ xx(i,:) = xopt((i-1)*Nmines+1:i*Nmines)';
+ Qt(i) = xopt(Nmines*Nyears+i);
+end
+years = string(1:Nyears);
+Mines = [];
+mprintf('Total profit %d',-fopt);
+for i = 1:Nmines
+ Mines = [Mines strcat(['Mine',string(i)])]
+end
+table = [['years' Mines];[years' string(xx)]];
+disp(table);
+disp(string(Qt),'Each year produced blended ore')
diff --git a/code/Symphonymat/factory.xls b/code/Symphonymat/factory.xls
new file mode 100644
index 0000000..04b174e
--- /dev/null
+++ b/code/Symphonymat/factory.xls
Binary files differ
diff --git a/code/fminbnd/Chichinadze.sce b/code/fminbnd/Chichinadze.sce
new file mode 100644
index 0000000..efd79f2
--- /dev/null
+++ b/code/fminbnd/Chichinadze.sce
@@ -0,0 +1,64 @@
+//Reference: Ernesto P. Adorio and U.P. Diliman,"MVF - Multivariate Test Functions
+//Library in C for Unconstrained Global Optimization",2005,
+//http://www.geocities.ws/eadorio/mvf.pdf, Last accessed, 11th June 2018
+//Example: The 2-dimensional function mvfChichinadze computes
+//f = x1^2 - 12*x1 + 11 + 10*cos(%pi/2*x1)) + 8*sin(5*%pi*x1)) - 1/sqrt(5)*exp(-(x2)-0.5)^2/2)
+//with domain −30 ≤ x0 ≤ 30, −10 ≤ x1 ≤ 10. The global minimum is 43.3159 at (5.90133, 0.5).
+//
+
+clc;
+//Objective Function
+function f = Chichinadze(x)
+ f = x(1)^2 - 12*x(1) + 11 + 10*cos(%pi/2*x(1)) + 8*sin(5*%pi*x(1)) - 1/sqrt(5)*exp(-(x(2)-0.5)^2/2)
+endfunction
+
+//Lower bound on the variables
+x1 = [-30 -10];
+//Upper bound on the variables
+x2 = [30 10];
+
+Maxit = 1500;
+CPU = 100;
+Tolx = 1e-6;
+
+mprintf('The termination criteria is as follows: Maximum Iterations = %d, Maximum CPU time = %d, Tolerance on solution = %f',Maxit,CPU,Tolx);
+
+//Options structure
+options=list("MaxIter",Maxit,"CpuTime", CPU,"TolX",Tolx)
+
+[xopt,fopt,exitflag,output,lambda]=fminbnd(Chichinadze,x1,x2,options)
+
+// Result representation
+clc;
+select exitflag
+case 0
+ disp("Optimal Solution Found")
+ disp(xopt',"The optimum solution obtained is")
+ disp(fopt,"The objective function value is")
+case 1
+ disp("Maximum Number of Iterations Exceeded. Output may not be optimal")
+ disp(xopt,"The solution obtained")
+ f = Chichinadze(xopt)
+ disp(f,"The objective function value is")
+case 2
+ disp("Maximum CPU Time exceeded. Output may not be optimal")
+ disp(xopt,"The solution obtained")
+ f = Chichinadze(xopt)
+ disp(f,"The objective function value is")
+case 3
+ disp("Stop at Tiny Step")
+ disp(xopt,"The solution obtained")
+ f = Chichinadze(xopt)
+ disp(f,"The objective function value is")
+case 4
+ disp("Solved To Acceptable Level")
+ disp(xopt,"The solution obtained")
+ f = Chichinadze(xopt)
+ disp(f,"The objective function value is")
+case 5
+ disp("Converged to a point of local infeasibility")
+ disp(xopt,"The solution obtained")
+ f = Chichinadze(xopt)
+ disp(f,"The objective function value is")
+end
+disp(output)
diff --git a/code/fminbnd/Cola.sce b/code/fminbnd/Cola.sce
new file mode 100644
index 0000000..203a3c4
--- /dev/null
+++ b/code/fminbnd/Cola.sce
@@ -0,0 +1,89 @@
+//Reference: Ernesto P. Adorio and U.P. Diliman,"MVF - Multivariate Test Functions
+//Library in C for Unconstrained Global Optimization",2005,
+//http://www.geocities.ws/eadorio/mvf.pdf, Last accessed, 11th June 2018
+//Example: The 17-dimensional function mvfCola computes indirectly the formula f(n,u)
+//by setting x0 = y0, x1 = u0, x(i) = u(2(i−2)), y(i) = u(2(i−2)+1)
+//f(n,u) = h(x,y) =sum((r(i,j) − d(i,j))^2) for all j<i
+//r(i,j) = [(x(i) − x(j) )^2 + (y(i) − y(j))^2]^1/2
+//This function has bounds 0 ≤ x0 ≤ 4 and −4 ≤ xi ≤ 4 for i = 1 . . . n − 1. It
+//has a global minimum of 11.7464.
+//u = [0.651906,1.30194,0.099242,-0.883791,-0.8796,0.204651,-3.28414,0.851188,-3.46245,2.53245,-0.895246,1.40992,-3.07367,1.96257,-2.97872,-0.807849,-1.68978];
+//
+
+clc;
+//Objective Function
+function summation = Cola(u)
+
+ d = [ 0 0 0 0 0 0 0 0 0 0
+ 1.27 0 0 0 0 0 0 0 0. 0
+ 1.69 1.43 0 0 0 0 0 0 0 0
+ 2.04 2.35 2.43 0 0 0 0 0 0 0
+ 3.09 3.18 3.26 2.85 0 0 0 0 0 0
+ 3.20 3.22 3.27 2.88 1.55 0 0 0 0 0
+ 2.86 2.56 2.58 2.59 3.12 3.06 0 0 0 0
+ 3.17 3.18 3.18 3.12 1.31 1.64 3.00 0 0 0
+ 3.21 3.18 3.18 3.17 1.70 1.36 2.95 1.32 0 0
+ 2.38 2.31 2.42 1.94 2.85 2.81 2.56 2.91 2.97 0];
+ summation = 0;
+
+ for i = 1:10
+ for j = 1:9
+ if j<i
+ x(1) = 0; y(1) = 0;y(2) = 0; x(2)= u(1);
+ if i>2
+ x(i) = u(2*(i-2));y(i) = u(2*(i-2)+1);
+ end
+ r(i,j) = ((x(i) - x(j))^2 + (y(i)-y(j))^2)^0.5;
+ summation = summation + (( r(i,j) - d(i,j) )^2)
+ end
+ end
+ end
+
+endfunction
+
+
+//Lower bound
+x1 = [0 -4*ones(1,16)];
+//Upper bound
+x2 = 4*ones(1,17);
+
+//Options structure
+options=list("MaxIter",[1500],"CpuTime", [100],"TolX",[1e-6])
+
+[xopt,fopt,exitflag,output,lambda]=fminbnd(Cola,x1,x2,options)
+
+// Result representation
+clc;
+select exitflag
+case 0
+ disp("Optimal Solution Found")
+ disp(xopt',"The optimum solution obtained is")
+ disp(fopt,"The objective function value is")
+case 1
+ disp("Maximum Number of Iterations Exceeded. Output may not be optimal")
+ disp(xopt,"The solution obtained")
+ f = Cola(xopt)
+ disp(f,"The objective function value is")
+case 2
+ disp("Maximum CPU Time exceeded. Output may not be optimal")
+ disp(xopt',"The solution obtained")
+ f = Cola(xopt)
+ disp(f,"The objective function value is")
+case 3
+ disp("Stop at Tiny Step")
+ disp(xopt',"The solution obtained")
+ f = Cola(xopt)
+ disp(f,"The objective function value is")
+case 4
+ disp("Solved To Acceptable Level")
+ disp(xopt',"The solution obtained")
+ f = Cola(xopt)
+ disp(f,"The objective function value is")
+case 5
+ disp("Converged to a point of local infeasibility")
+ disp(xopt',"The solution obtained")
+ f = Cola(xopt)
+ disp(f,"The objective function value is")
+end
+disp(output)
+
diff --git a/code/fmincon/Tankprob.sci b/code/fmincon/Tankprob.sci
new file mode 100644
index 0000000..500090e
--- /dev/null
+++ b/code/fmincon/Tankprob.sci
@@ -0,0 +1,95 @@
+//Design a circular tank, closed at both ends, with a volume of 200 m3.The cost is proportional to the surface area of material, which is priced
+//at $400/m2. The tank is contained within a shed with a sloping roof,thus the height of the tank h is limited by h ≤ 12 − d/2, where d is
+//the tank diameter. Formulate the minimum cost problem and solve the design problem.
+//Ref:Diwekar, Urmila,Introduction to Applied Optimization, Introduction to Applied Optimization, Editor:Ding-Zhu Du, Springer Optimization and Its Applications Springer Optimization and Its Applications, VOL 22, Chapter 3
+
+// 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;
+
+function cost = Tankprob(x)
+ r = x(1);h = x(2);
+ cost = (2*%pi*r*h+2*%pi*r^2)*400;
+endfunction
+
+function [ceq,c] = Nonlinearcon(x)
+ r = x(1);h = x(2);
+ c = [];
+ ceq = 200-%pi*r^2*h;
+endfunction
+
+function y=Gradobj(x)
+ y= [800*%pi*x(2) + 1600*%pi*x(1),800*%pi*x(1)];
+endfunction
+
+mprintf('\nDesign a circular tank, closed at both ends, with a volume of 200 m3 with minimum cost.\n The tank is contained within a shed with a sloping roof,thus the height of the tank h is limited')
+mprintf('\nCost of material is: %f',400);
+mprintf('The design variables are radius and height of the tank')
+
+A = [1 1];b = 12;
+
+x0 = input('Enter initial guess as vector:');
+if (sum(x0<=0) | (length(x0)~=2))
+ x0 = [1 2];
+ mprintf('Incorrect initial guess...\n changing initial guess to r = %d and h = %d',x0(1),x0(2));
+end
+
+
+lb = [0 0];
+input('press enter to continue')
+options=list("MaxIter",1000,"GradObj", Gradobj);
+
+[xopt,fopt,exitflag,output1] = fmincon(Tankprob,x0,A,b,[],[],lb,[],Nonlinearcon,options);
+
+[xopt1,fopt1,exitflag1,output2] = fmincon(Tankprob,x0,A,b,[],[],lb,[],Nonlinearcon);
+
+clc
+select exitflag
+case 0
+ mprintf('Optimal Solution Found');
+ mprintf('\nThe optimal radius and height of the tank are %f m and %f m',xopt(1),xopt(2));
+ mprintf('\nThe volume of the tank is %f m^3',%pi*xopt(1)^2*xopt(2));
+ mprintf('\nThe total surface area and cost of the tank are %f m^2 and %f $',fopt/400,fopt)
+ mprintf('\nTime taken to solve the problem with gradient information is %f s and without gradient information is %f s',output1.Cpu_Time,output2.Cpu_Time);
+case 1
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal');
+ input('press enter to view results');
+ printf('\nThe optimal radius and height of the tank are %f m and %f m',xopt(1),xopt(2));
+ mprintf('\nThe volume of the tank is %f m^3',%pi*xopt(1)^2*xopt(2));
+ mprintf('\nThe total surface area and cost of the tank are %f m^2 and %f $',fopt/400,fopt)
+case 2
+ mprintf('Maximum amount of CPU Time exceeded. Output may not be optimal.');
+ input('press enter to view results');
+ printf('\nThe optimal radius and height of the tank are %f m and %f m',xopt(1),xopt(2));
+ mprintf('\nThe volume of the tank is %f m^3',%pi*xopt(1)^2*xopt(2));
+ mprintf('\nThe total surface area and cost of the tank are %f m^2 and %f $',fopt/400,fopt)
+case 3
+ mprintf('Stop at Tiny Step');
+ input('press enter to view results');
+ printf('\nThe optimal radius and height of the tank are %f m and %f m',xopt(1),xopt(2));
+ mprintf('\nThe volume of the tank is %f m^3',%pi*xopt(1)^2*xopt(2));
+ mprintf('\nThe total surface area and cost of the tank are %f m^2 and %f $',fopt/400,fopt)
+case 4
+ mprintf('Solved To Acceptable Level');
+ input('press enter to view results');
+ printf('\nThe optimal radius and height of the tank are %f m and %f m',xopt(1),xopt(2));
+ mprintf('\nThe volume of the tank is %f m^3',%pi*xopt(1)^2*xopt(2));
+ mprintf('\nThe total surface area and cost of the tank are %f m^2 and %f $',fopt/400,fopt)
+case 5
+ mprintf('Converged to a point of local infeasibility.');
+ input('press enter to view results');
+ printf('\nThe optimal radius and height of the tank are %f m and %f m',xopt(1),xopt(2));
+ mprintf('\nThe volume of the tank is %f m^3',%pi*xopt(1)^2*xopt(2));
+ mprintf('\nThe total surface area and cost of the tank are %f m^2 and %f $',fopt/400,fopt)
+end
+
+
diff --git a/code/fmincon/spring.sci b/code/fmincon/spring.sci
new file mode 100644
index 0000000..0409dce
--- /dev/null
+++ b/code/fmincon/spring.sci
@@ -0,0 +1,120 @@
+//The problem consists of minimizing the weight of a tension/compression spring subject to constraints on
+//minimum deflection, shear stress, surge frequency, limits on outside diameter and on design variables. The design
+//variables are the mean coil diameter D, the wire diameter d and the number of active coils N. The problem can be
+//expressed as follows:
+//Min f(x) = ( N + 2)*D*d^2
+//ST
+//g1(x) = 1-D^3*N/(71785*d^4) <= 0;
+//g2(x) = 1-140.45*d/(D^2*N) <= 0;
+//g3(x) = (4*D^2-d*D)/(12566*(d^3*D-d^4)) + 1/(5108*d^2)-1 <= 0;
+//g4(x) = (D+d)/1.5 - 1 <=0;
+//The following ranges for the variables were used
+//0.05 <= d <= 2.0, 0.25<= D <=1.3, 2.0<=N<=15.0
+//
+//Ref:Arora, J. S.. Introduction IO optimum design New York McGraw-Hill, 1989
+//======================================================================
+// 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;
+
+function y = spring(x)
+ y = (x(3)+2)*x(1)^2*x(2);
+endfunction
+
+function [C,Ceq] = Nonlinconstraint(x)
+ C(1) = 1-x(2)^3*x(3)/(71785*x(1)^4);
+ C(2) = 1-140.45*x(1)/(x(2)^2*x(3));
+ C(3) = (4*x(2)^2-x(1)*x(2))/(12566*(x(1)^3*x(2)-x(1)^4)) + 1/(5108*x(1)^2)-1;
+ Ceq = [];C = C';
+endfunction
+
+function f = fGrad(x)
+ f = [2*x(1)*x(2)*(x(3)+2) x(1)^2*(x(3)+2) x(1)^2*x(2)];
+endfunction
+
+mprintf('minimizing the weight of a tension/compression spring ');
+mprintf('\n\n Design variables include mean coil diameter:D, the wire diameter:d and the number of active coils:N ');
+disp('The objective is :( N + 2)*D*d^2');
+mprintf('\nNon-linear constraints are on minimum deflection, shear stress, surge frequency,');
+mprintf('\nlinear constraint is on outside diameter');
+
+
+A = [1 1 0];
+b = 1.5;
+mprintf('Bounds on the variable are as follows');
+
+lb = [0.05 0.25 2];
+ub = [2 1.3 15];
+table = [['bounds', 'lb','ub'];[['d';'D';'N'],string(lb'),string(ub')]];
+disp(table);
+x0 = lb + rand(1)*(ub-lb);
+disp(x0,'Initial Guess is :');
+
+input('press enter to continue');
+
+mprintf('Scilab is solving your problem');
+options=list("MaxIter", [1500], "CpuTime", [500], "GradObj", fGrad);
+
+[xopt,fval,exitflag,output] =fmincon(spring, x0,A,b,[],[],lb,ub,Nonlinconstraint,options);
+
+clc
+select exitflag
+case 0
+ mprintf('Optimal Solution Found');
+ mprintf('\n\nThe optimal design parameters are d = %f ,D = %f and N = %f m',xopt(1),xopt(2),xopt(3));
+ mprintf('\n\nThe optimal objective is %f ',fval);
+ mprintf('\n\nIterations completed %d',output.Iterations);
+ mprintf('\nTotal CPU time %f',output.Cpu_Time);
+ mprintf('\nTotal Functional Evaluations %d',output.Objective_Evaluation);
+case 1
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal');
+ input('press enter to view results');
+ mprintf('\nThe optimal design parameters are d = %f ,D = %f and N = %f m',xopt(1),xopt(2),xopt(3));
+ mprintf('\nThe optimal objective is %f ',fval);
+ mprintf('\nIterations completed %d',output.Iterations);
+ mprintf('\nTotal CPU time %f',output.Cpu_Time);
+ mprintf('\nTotal Functional Evaluations %d',output.Objective_Evaluation);
+case 2
+ mprintf('Maximum amount of CPU Time exceeded. Output may not be optimal.');
+ input('press enter to view results');
+ mprintf('\nThe optimal design parameters are d = %f ,D = %f and N = %f m',xopt(1),xopt(2),xopt(3));
+ mprintf('\nThe optimal objective is %f ',fval);
+ mprintf('\nIterations completed %d',output.Iterations);
+ mprintf('\nTotal CPU time %f',output.Cpu_Time);
+ mprintf('\nTotal Functional Evaluations %d',output.Objective_Evaluation);
+case 3
+ mprintf('Stop at Tiny Step');
+ input('press enter to view results');
+ mprintf('\nThe optimal design parameters are d = %f ,D = %f and N = %f m',xopt(1),xopt(2),xopt(3));
+ mprintf('\nThe optimal objective is %f ',fval);
+ mprintf('\nIterations completed %d',output.Iterations);
+ mprintf('\nTotal CPU time %f',output.Cpu_Time);
+ mprintf('\nTotal Functional Evaluations %d',output.Objective_Evaluation);
+case 4
+ mprintf('Solved To Acceptable Level');
+ input('press enter to view results');
+ mprintf('\nThe optimal design parameters are d = %f ,D = %f and N = %f m',xopt(1),xopt(2),xopt(3));
+ mprintf('\nThe optimal objective is %f ',fval);
+ mprintf('\nIterations completed %d',output.Iterations);
+ mprintf('\nTotal CPU time %f',output.Cpu_Time);
+ mprintf('\nTotal Functional Evaluations %d',output.Objective_Evaluation);
+case 5
+ mprintf('Converged to a point of local infeasibility.');
+ input('press enter to view results');
+ mprintf('\nThe optimal design parameters are d = %f ,D = %f and N = %f m',xopt(1),xopt(2),xopt(3));
+ mprintf('\nThe optimal objective is %f ',fval);
+ mprintf('\nIterations completed %d',output.Iterations);
+ mprintf('\nTotal CPU time %f',output.Cpu_Time);
+ mprintf('\nTotal Functional Evaluations %d',output.Objective_Evaluation);
+end
+
+
+
diff --git a/code/fminsearch/Brownsfunc.sci b/code/fminsearch/Brownsfunc.sci
new file mode 100644
index 0000000..e64a650
--- /dev/null
+++ b/code/fminsearch/Brownsfunc.sci
@@ -0,0 +1,73 @@
+// This is an example for unconstraint nonlinear problems.
+//Ref:J. J. More, B. S. Garbow, and K. E. Hillstrom, Testing unconstrained optimization software, ACM Transactions on Mathematical Software, Vol. 7, No. 1, pp. 17–41, 1981.
+//Example:
+//f(x1,x2) = (x1 - 10^6)^2 + (x2 - 2*10^-6)^2 + (x1*x2 - 2)^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;close
+
+function y = Brownsfunc(x)
+ y = (x(1)-1d6)^2 + (x(2)-2*1D-6)^2 + (x(1)*x(2)-2)^2
+endfunction
+
+function stop=outfun(x, optimValues, state)
+ subplot(1,2,1)
+ plot(optimValues.funccount,optimValues.fval,'r.');
+ xlabel('function count');ylabel('Objective value')
+
+ subplot(1,2,2)
+ plot(optimValues.funccount,x(1),'r.');
+ plot(optimValues.funccount,x(2),'b.');
+ legend(['X1','X2'])
+ set(gca(),"auto_clear","off")
+ xlabel('function count');ylabel('variable values')
+
+ stop = %f
+endfunction
+
+X0 = [1 1];
+MFes = 500;
+Miter = 200;
+TF = 1D-6;
+TX = 1D-6;
+mprintf('The following settings are used\n Maximum iterations %d \n maximum functional exaluations %d\n Function tolerance %s \n variable tolerance %s ',Miter,MFes,string(TF),string(TX));
+mprintf('\nThe initial guess is x1 = %f and x2 = %f',X0(1),X0(2))
+input('Press enter to proceed ')
+clc;
+mprintf('Scilab is solving the problem...')
+
+options = optimset ("MaxFunEvals",MFes,"MaxIter",Miter,"TolFun",TF,"TolX",TX, "OutputFcn" , outfun);
+
+[x,fval,exitflag,output] = fminsearch(Brownsfunc,X0,options)
+
+clc
+select exitflag
+case -1
+ disp(output.algorithm, 'Algorithm used')
+ mprintf('\n The maximum number of iterations has been reached \n')
+ mprintf('\n The number of iterations %d ',output.iterations)
+ mprintf('\n The number of function evaluations %d',output.funcCount)
+case 0
+ disp(output.algorithm, 'Algorithm used ')
+ mprintf('\n The maximum number of function evaluations has been reached \n')
+ mprintf('\n The number of function evaluations %d',output.funcCount)
+ mprintf('\n The number of iterations %d ',output.iterations)
+
+case 1
+ disp(output.algorithm, 'Algorithm used ')
+ mprintf('\n The tolerance on the simplex size and function value delta has been reached\n')
+ mprintf('\n The number of function evaluations %d',output.funcCount)
+ mprintf('\n The number of iterations %d ',output.iterations)
+end
+
+disp(x,"The optimal solution is")
+mprintf("\n The optimum value of the function is %s",string(fval))
diff --git a/code/fminsearch/FletcherPowell.sce b/code/fminsearch/FletcherPowell.sce
new file mode 100644
index 0000000..b3b4971
--- /dev/null
+++ b/code/fminsearch/FletcherPowell.sce
@@ -0,0 +1,65 @@
+// This is an example for unconstraint nonlinear problems.
+// Ref:R.fletcher and M.J.D Powell, A Rapidly Convergent Descent Method for Minimization Algorithms, Computer journal, Vol. 6, pp. 163-168, 1963
+//Example:
+//f(x1,x2,x3) = 100*((x3 - 10*theta(x1,x2))^2 + (sqrt(x1^2 + x1^2) - 1)^2) + x3^2
+//theta(x1,x2) = (atan(x(2)/x(1)))/(2*%pi) if x(1)>0
+// = %pi + atan(x(2)/x(1)) if x(1)<0
+//======================================================================
+// 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;
+
+clc;clear;close
+
+function y = FletcherPowell(x)
+ if (x(1)>0)
+ theta_x1x2 = (atan(x(2)/x(1)))/(2*%pi);
+ elseif (x(1)<0)
+ theta_x1x2 = %pi + atan(x(2)/x(1));
+ end
+ y = 100*( (x(3) - 10*theta_x1x2 ).^2 + (sqrt(x(1)^2 + x(2)^2) - 1)^2) + x(3)^2;
+endfunction
+
+X0 = [-1 0 0];
+MFes = 500;
+Miter = 200;
+TF = 1D-10;
+TX = 1D-10;
+mprintf('The following settings are used\n Maximum iterations %d \n maximum functional exaluations %d\n Function tolerance %s \n variable tolerance %s ',Miter,MFes,string(TF),string(TX));
+input('Press enter to proceed ')
+clc;
+mprintf('Scilab is solving the problem...')
+
+options = optimset ("MaxFunEvals",MFes,"MaxIter",Miter,"PlotFcns",optimplotfval,"TolFun",TF,"TolX",TX);
+
+[x,fval,exitflag,output] = fminsearch(FletcherPowell,X0,options)
+clc
+select exitflag
+case -1
+ disp(output.algorithm, 'Algorithm used')
+ mprintf('\n The maximum number of iterations has been reached \n')
+ mprintf('\n The number of iterations %d ',output.iterations)
+ mprintf('\n The number of function evaluations %d',output.funcCount)
+case 0
+ disp(output.algorithm, 'Algorithm used ')
+ mprintf('\n The maximum number of function evaluations has been reached \n')
+ mprintf('\n The number of function evaluations %d',output.funcCount)
+ mprintf('\n The number of iterations %d ',output.iterations)
+
+case 1
+ disp(output.algorithm, 'Algorithm used ')
+ mprintf('\n The tolerance on the simplex size and function value delta has been reached\n')
+ mprintf('\n The number of function evaluations %d',output.funcCount)
+ mprintf('\n The number of iterations %d ',output.iterations)
+end
+
+disp(x,"The optimal solution is")
+mprintf("\n The optimum value of the function is %s",string(fval))
diff --git a/code/fminunc/Freudenstein and Roth.sci b/code/fminunc/Freudenstein and Roth.sci
new file mode 100644
index 0000000..a720c97
--- /dev/null
+++ b/code/fminunc/Freudenstein and Roth.sci
@@ -0,0 +1,75 @@
+// This is an example of unconstrained optimization problem using fminunc function. The problem has global optima at x* = [5 4] and f* = 0;
+//An alternate optima to the problem is x = [11.41.. -0.8968..] and f = 48.9842..
+//Ref:F. Freudenstein and B. Roth, Numerical solution of systems of nonlinear equations, Journal of ACM , Vol. 10, No. 4, pp. 550–556, 1963.
+//Ref:S.S. Rao, “Engineering optimization: Theory and Practice”, John Wiley & Sons Inc., New York (NY), 3rd edition edition, 1996. Chapter 6
+//======================================================================
+// 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;
+
+function y = Freud_Roth(x)
+ y = (-13 + x(1) + ( (5-x(2)) *x(2) -2 )*x(2) )^2 + ( -29 + x(1) + (( x(2)+1 )*x(2) -14 )*x(2) )^2;
+endfunction
+
+mprintf('Freudenstein and Roth function: This function has 2 variables')
+//Initial point
+x0 = input('Enter initial point in vector form ' )
+if length(x0)~=2
+ mprintf('Incorrect initial point. Taking initial point as [0.5 -2]');
+ x0 = [0.5 -2];
+end
+
+//Gradient of objective function using numderivative function
+function y=Grad(x)
+y = numderivative(Freud_Roth,x);
+endfunction
+
+//Hessian of Objective Function using numderivative function
+function y=Hess(x)
+[G,H] = numderivative(Freud_Roth,x);
+Nvar = length(x);
+for i = 1:Nvar
+ y(i,:) = H((i-1)*Nvar+1:i*Nvar);
+end
+endfunction
+
+//Options structure
+options=list("MaxIter", [1500], "CpuTime", [500], "gradobj", Grad, "hessian", Hess);
+//Calling Ipopt
+[xopt,fopt,exitflag,output,gradient,hessian]=fminunc(Freud_Roth,x0,options)
+clc
+select exitflag
+case 0
+ disp("Optimal Solution Found")
+ disp(xopt', "The optimum solution obtained is")
+ disp(fopt, "The optimum objective function value is")
+case 1
+ disp("Maximum Number of Iterations Exceeded. Output may not be optimal.")
+ disp(xopt', "The solution obtained is")
+ disp(fopt, "The objective function value is")
+case 2
+ disp("Maximum CPU Time exceeded. Output may not be optimal.")
+ disp(xopt', "The solution obtained is")
+ disp(fopt, "The objective function value is")
+case 3
+ disp("Stop at Tiny Step.")
+ disp(xopt', "The solution obtained is")
+ disp(fopt, "The objective function value is")
+case 4
+ disp("Solved To Acceptable Level.")
+ disp(xopt', "The solution obtained is")
+ disp(fopt, "The objective function value is")
+case 5
+ disp("Converged to a point of local infeasibility.")
+end
+disp(output)
+
+
diff --git a/code/fsolve/Redlich-Kwong.sce b/code/fsolve/Redlich-Kwong.sce
new file mode 100644
index 0000000..c2fc326
--- /dev/null
+++ b/code/fsolve/Redlich-Kwong.sce
@@ -0,0 +1,69 @@
+//Ref:Steven C. Chapra. 2006. Applied Numerical Methods with MATLAB for Engineers and Scientists. McGraw-Hill Science/Engineering/Math,Chapter 6
+//Example:
+//The Redlich-Kwong equation of state is given by
+//p = ((R*T)/(v-b) - a/(v*(v+b)*sqrt(T)))
+//where R = the universal gas constant [= 0.518 kJ/(kg K)], T = absolute temperature (K), p = absolute pressure (kPa),and v = the volume of a kg of gas (m3/kg). The parameters a and b are calculated by
+// a = 0.427*(R^2*Tc^2.5)/pc; b = 0.0866*R*(Tc/pc);
+//where pc = 4600 kPa and Tc = 191 K. As a chemical engineer, you are asked to determine the amount of methane fuel that can be held in a 3 m3 tank at a temperature of −40 ◦C with a pressure of 65,000 kPa. Use a root-locating method of your choice to calculate v and then determine the mass of methane contained in the tank.
+//Note: The initial guess has to be assumed by the user. An improper initial guess will result in deviation from the root.
+//======================================================================
+// 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;
+
+function y = Redlich_Kwong(v)
+
+ R = 0.518;T = 273+(-40);p = 65000;Tc = 191;pc = 4600;
+ a = 0.427*(R^2*Tc^2.5)/pc; b = 0.0866*R*(Tc/pc);
+
+ y = p - ( (R*T)/(v-b) - a/(v*(v+b)*sqrt(T)) );
+
+endfunction
+
+function dy = dv_Redlich_Kwong(x)
+
+ dy = numderivative(Redlich_Kwong, x);
+
+endfunction
+
+// Set of initial values to check
+x0 = [1 0.01 -1.0 -0.01] ;
+
+tol = 1D-10;
+
+for i =1:length(x0)
+
+ disp(x0(i),'Initial guess value is')
+
+ [x ,v ,info]=fsolve(x0(i),Redlich_Kwong ,dv_Redlich_Kwong ,tol)
+
+ select info
+ case 0
+ mprintf('\n improper input parameters\n');
+ case 1
+ mprintf('\n algorithm estimates that the relative error between x and the solution is at most tol\n');
+ case 2
+ mprintf('\n number of calls to fcn reached\n');
+ case 3
+ mprintf('\n tol is too small. No further improvement in the approximate solution x is possible\n');
+ else
+ mprintf('\n iteration is not making good progress\n');
+ end
+
+ mprintf('\n volume of metheane is %f m3/kg and Mass of methane is %f kg\n',x,3/x);
+
+ if i<length(x0)
+ input('press enter to check the next initial guess value')
+ clc
+ end
+
+end
diff --git a/code/fsolve/SLEs.sce b/code/fsolve/SLEs.sce
new file mode 100644
index 0000000..3daa1ee
--- /dev/null
+++ b/code/fsolve/SLEs.sce
@@ -0,0 +1,50 @@
+//This problem solves sustem of nonlinear equations using fsolve
+//Reference: A. Golbabai, M. Javidi,Newton-like iterative methods for solving system of non-linear equations,Applied Mathematics and Computation,Volume 192, Issue 2,2007,Pages 546-551,ISSN 0096-3003,https://doi.org/10.1016/j.amc.2007.03.035.(http://www.sciencedirect.com/science/article/pii/S0096300307003578)
+//======================================================================
+// 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
+//Function representing system of nonlinear equaltions
+function f = SLEs(p)
+
+ x = p(1);y = p(2);
+ f(1) = (x+10)/x^4 + x^2 + 3*x*y + y^2-16;
+ f(2) = 1/(x+y)^7 + x*y - 1 - 1/2^7;
+
+endfunction
+
+//Tolarence of solution
+tol = 1D-15;
+
+//Initial guess
+x0 = [10 10];
+
+disp(x0,'Initial guess value is');
+
+//Obtaining solution using fsolve
+[x ,v ,info] = fsolve(x0,SLEs,tol)
+clc
+select info
+case 0
+ mprintf('\n improper input parameters\n');
+case 1
+ mprintf('\n algorithm estimates that the relative error between x and the solution is at most tol\n');
+case 2
+ mprintf('\n number of calls to fcn reached\n');
+case 3
+ mprintf('\n tol is too small. No further improvement in the approximate solution x is possible\n');
+else
+ mprintf('\n iteration is not making good progress\n');
+end
+disp(x,'The solution is ')
+disp(v,'the function value at solution')
+
+
diff --git a/code/intfmincon/ProcessSel.sce b/code/intfmincon/ProcessSel.sce
new file mode 100644
index 0000000..d2efefb
--- /dev/null
+++ b/code/intfmincon/ProcessSel.sce
@@ -0,0 +1,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
+
+
+
diff --git a/code/intqpipopt/SpringCart.sce b/code/intqpipopt/SpringCart.sce
new file mode 100644
index 0000000..c676173
--- /dev/null
+++ b/code/intqpipopt/SpringCart.sce
@@ -0,0 +1,68 @@
+//Three carts interconnected by springs, are subjected to the loads P1, P2, and P3. The displacementas of the carts can be found by minimizing the potential energy of the system.
+//
+//f = 0.5X'[K]X - X'P
+//
+//where [K] = [k1+k4+k5 -k4 - k5
+// -k4 k(2)+k(4)+k(6) -k6
+// -k5 -k(6) k3+k5+k6+k7+k8 ]
+//The following data is given: k1 = 5000 N/m , k2 = 1500 N/m,
+//k3 = 2000 N/m, k4 = 1000 N/m, k5 = 2500 N/m, k6 = 500 N/m, k7 = 3000 N/m, k8 = 3500 N/m
+//P1 = 1000 N/m, P2 = 2000 N/m and P3 = 3000 N/m.
+//The optimum displacement is x = [0.3241 0.8360 0.3677];
+//Reference: Rao, S. S. (2009). Engineering Optimization: Theory and Practice: Fourth Edition. John Wiley and Sons. Chapter 6.DOI: 10.1002/9780470549124.
+
+//======================================================================
+// 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
+// Given data
+k = [5000 1500 2000 1000 2500 500 3000 3500];
+p = [1000 2000 3000];
+
+Kval = [k(1)+k(4)+k(5) -k(4) -k(5)
+ -k(4) k(2)+k(4)+k(6) -k(6)
+ -k(5) -k(6) k(3)+k(5)+k(6)+k(7)+k(8) ]
+
+
+//Initial Guess
+x0 = [1 1 1];
+
+//Unconstraint,unbounded and continuous problem. Hence the values of intcon,bounds and constraints are giben empty matrices
+//Setting Iteration and cpu time in options
+options = list("MaxIter", 100, "CpuTime", 100);
+
+//calling the optimization solver
+[xopt,fopt,exitflag,output] = intqpipopt(Kval,-p,[],[],[],[],[],[],[],x0,options)
+
+clc
+//Result display
+select exitflag
+case 0
+disp("Optimal Solution Found")
+printf('\n The displacement of the three springs are \n x1 = %f\n x2 = %f\n x3 = %f',xopt(1),xopt(2),xopt(3));
+printf('\nThe potential energy of the system is %f',fopt)
+printf('\nThe constraint violation is %f',output.ConstrViolation)
+case 1
+ disp("InFeasible Solution")
+ printf('\nThe constraint violation is %f',output.ConstrViolation)
+case 2
+ disp("Objective Function is Continuous Unbounded.")
+ printf('\nThe constraint violation is %f',output.ConstrViolation)
+case 3
+ disp("Limit Exceeded.")
+ printf('\nThe constraint violation is %f',output.ConstrViolation)
+case 4
+ disp("User Interrupt.")
+case 5
+ disp("MINLP Error.")
+end
+
diff --git a/code/linprog/Blending_problem.sce b/code/linprog/Blending_problem.sce
new file mode 100644
index 0000000..24b80ad
--- /dev/null
+++ b/code/linprog/Blending_problem.sce
@@ -0,0 +1,135 @@
+// A practcal problem of blending of different products to acheive designed performance levels.
+// This is an example from the class of allocation problems.
+// Ref:H.A. Eiselt and C.-L.SAndblom,"Linear Programming and its Applications",Springer-Verlag Berlin Heidelberg 2007,chapter 2.7
+//Example:
+//A firm faces the problem of blending three raw materials into two final
+//products. The required numerical information is provided in Table.
+//
+//-----------------------------------------------------------------------
+// Raw Products Amount available Unit Cost
+//Materials 1 2 (in tons) ($)
+//-----------------------------------------------------------------------
+// 1 [0.4;0.6][0.5;0.6] 2000 1
+// 2 [0.1;0.2][0.1;0.4] 1000 1.5
+// 3 [0.2;0.5][0.2;0.3] 500 3
+//-----------------------------------------------------------------------
+//quantity required 600 700
+//(in tons)
+//-----------------------------------------------------------------------
+//unit selling price 10 8
+//($)
+//-----------------------------------------------------------------------
+//we assume that raw materials blend linearly, meaning that taking,
+//say, α units of raw material A and β units of raw material B, then
+//the resulting blend C has features that are proportional to the
+//quantities of A and B that C is made of. As an example, take 3 gallons
+//of 80º water and 2 gallons of 100º water,then the result would be 5
+//gallons of water, whose temperature is [3(80) + 2(100)]/5 = 88º.
+//The parameters max and min indicate the smallest and largest proportion
+//of rawmaterial that is allowed in the final product.Determine the
+//optimal transportation plan to maximize profit.
+
+// 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;
+Nprod = 2;
+Nraw = 3;
+demand = [600 700];
+cost = [10 8];
+con(1,:,:) = [0.4 0.6;0.1 0.2;0.2 0.5];
+con(2,:,:) = [0.5 0.6;0.1 0.4;;0.2 0.3];
+raw = [2000 1000 500];
+RawC = [1 1.5 3];
+mprintf('Received Data')
+mprintf('Number of products %d',Nprod);
+mprintf('number of raw materials %d',Nraw);
+disp(demand,'Demand of products');
+disp(cost,'Unit cost of products');
+
+product = [];Min_Max = [];
+for i = 1:Nprod
+ product = [product,['Min-Max',string(i)]];
+ for j = 1:Nraw
+ MM(j,:) = con(i,j,:);
+ end
+ Min_Max = [Min_Max string(MM)];
+end
+Row1 = ['Raw Material',product,'Amount_avail','unit_Cost']
+RawMaterial = string([1:Nraw]');
+Amount_avail = string(raw(:));
+unit_C = string(RawC(:));
+Table = [Row1;[RawMaterial Min_Max Amount_avail unit_C]];
+disp(Table);
+input('Press Enter to proceed ')
+A1 = zeros(Nraw,Nraw*Nprod);
+for i = 1:Nraw
+ A1(i,i:Nraw:Nraw*Nprod)=1;
+ b1(i) = raw(i);
+end
+A2 = zeros(Nprod*Nraw*2,Nprod*Nraw);b2 = zeros(Nprod*Nraw*2,1);
+A3 = zeros(Nprod,Nprod*Nraw);
+for i = 1:Nprod
+ Areq = ones(2*Nraw,Nraw);
+ Areq(Nraw+1:2*Nraw,:) = -1;
+ for j =1:Nraw
+ Areq(j,:) = con(i,j,1)* Areq(j,:);
+ Areq(j,j) = Areq(j,j)-1;
+ Areq(Nraw+j,:) = Areq(Nraw+j,:)*con(i,j,2);
+ Areq(Nraw+j,j) = 1 + Areq(Nraw+j,j);
+ end
+ A2((i-1)*2*Nraw+1:i*2*Nraw,(i-1)*Nraw+1:i*Nraw) = Areq;
+ A3(i,(i-1)*Nraw+1:i*Nraw)= ones(1,Nraw);
+ b3(i) = demand(i);
+end
+
+A = [A1;A2];b = [b1;b2];
+Aeq = A3;beq = b3;
+C = [];
+for i = 1:Nprod
+ C1 = -cost(i)*ones(1,Nraw) + RawC;
+ C = [C ;C1'];
+end
+
+lb = zeros(1,Nprod*Nraw);
+intcon = 1:Nprod*Nraw;
+[xopt,fopt,exitflag,output,lambda] = linprog(C,A,b,Aeq,beq,lb,[]);
+clc
+select exitflag
+case 0 then
+ mprintf('Optimal Solution Found')
+ input('Press enter to view results')
+ //Display results
+ mprintf('Total profit is %d\n',-fopt);
+
+ A1 = [];product = [];
+ for i=1:Nprod
+ Raw_M(:,i) = xopt((i-1)*Nraw+1:i*Nraw) ;
+ product(i) = string(sum(Raw_M(:,i)));
+ A1 = [A1 strcat(['Product' string(i)])];
+ end
+
+ table = [['RawMaterial', A1];[RawMaterial string(Raw_M)];['Total Product' product']];
+ disp(table);
+case 1 then
+ mprintf('Primal Infeasible')
+case 2 then
+ mprintf('Dual Infeasible')
+case 3
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal')
+case 4
+ mprintf('Solution Abandoned')
+case 5
+ mprintf('Primal objective limit reached')
+else
+ mprintf('Dual objective limit reached')
+end
+
+
diff --git a/code/linprog/Employee Scheduling.sce b/code/linprog/Employee Scheduling.sce
new file mode 100644
index 0000000..5fdce02
--- /dev/null
+++ b/code/linprog/Employee Scheduling.sce
@@ -0,0 +1,83 @@
+//Employee Scheduling
+//Macrosoft has a 24-hour-a-day, 7-days-a-week toll free hotline that is being set up to answer questions regarding a new product. The following table summarizes the number of full-time equivalent employees (FTEs) that mustbe on duty in each time block.
+//
+//Interval Time FTEs
+// 1 0-4 15
+// 2 4-8 10
+// 3 8-12 40
+// 4 12-16 70
+// 5 16-20 40
+// 6 20-0 35
+//Macrosoft may hire both full-time and part-time employees. The former work 8-hour shifts and the latter work 4-hour shifts; their respective hourly wages are $15.20 and $12.95. Employees may start work only at the beginning of 1 of the 6 intervals.Part-time employees can only answer 5 calls in the time a full-time employee can answer 6 calls. (i.e., a part-time employee is only 5/6 of a full-time employee.) At least two-thirds of theemployees working at any one time must be full-time employees.Formulate an LP to determine how to staff the hotline at minimum cost.
+
+// 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;
+
+Interval = 1:6;
+Time = 0:4:20;
+FTEs = [15 10 40 70 40 35];
+intr = [];
+Ninterval = length(Interval);
+for i = 1:Ninterval-1
+ intr = [intr;strcat([string(Time(i)),'-',string(Time(i+1))])];
+end
+
+intr = [intr;strcat([string(Time(i)),'-',string(Time(1))])];
+
+table = [['Interval','Time','Full-time Emp Req'];[string(Interval'),intr,string(FTEs')]];
+mprintf('Data received');
+disp(table);
+input('press enter to proceed');
+
+A1 = zeros(Ninterval,2*Ninterval);b1 = zeros(Ninterval,1);
+A2 = zeros(Ninterval,2*Ninterval);b2 = zeros(Ninterval,1);
+A1(1,[1 Ninterval Ninterval+1]) = [-1 -1 -5/6]; b1(1) = -15;
+A2(1,[1 Ninterval Ninterval+1]) = [-1/3 -1/3 2/3];
+for i = 2:Ninterval
+ A1(i,i-1:i) = -1;
+ A1(i,Ninterval+i) = -5/6;
+ b1(i) = -FTEs(i);
+
+ A2(i,i-1:i)=-1/3;
+ A2(+i,Ninterval+i)=2/3;
+end
+A= [A1;A2];b = [b1;b2];
+Cost = [15.20*8*ones(1,Ninterval) 12.95*4*ones(1,Ninterval)]';
+lb = zeros(1,2*Ninterval);ub = [];
+[xopt,fopt,exitflag,output,lambda] = linprog(Cost,A,b,[],[],lb,ub)
+
+clc
+select exitflag
+case 0 then
+ mprintf('Optimal Solution Found')
+ input('Press enter to view results')
+ //Display results
+ mprintf('Total Cost is %d\n',-fopt);
+ table(:,4:5) = [['Full time Emp', 'Part-Time Emp'];[string(xopt(1:6)) string(xopt(7:12))]];
+ disp(table);
+case 1 then
+ mprintf('Primal Infeasible')
+case 2 then
+ mprintf('Dual Infeasible')
+case 3
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal')
+case 4
+ mprintf('Solution Abandoned')
+case 5
+ mprintf('Primal objective limit reached')
+else
+ mprintf('Dual objective limit reached')
+end
+
+
+
diff --git a/code/linprog/FuelOil.sce b/code/linprog/FuelOil.sce
new file mode 100644
index 0000000..305cd7a
--- /dev/null
+++ b/code/linprog/FuelOil.sce
@@ -0,0 +1,86 @@
+// This example shows the use of equality constraint and corrosponding unrestricted dual variables
+//Example:
+//A refinery produces, on average, 1000 gallon/hour of virgin pitch in its crude distillation operation. This pitch may be blended with flux stock to make commercial fuel oil, or it can be sent in whole or in part to a visbreaker unit. The visbreaker produces an 80 percent yield of tar that can also be blended with flux stock to make commercial fuel oil. The visbreaking operation is economically break-even if the pitch and the tar are given no value, that is, the value of the overhead product equals the cost of the operation. The commercial fuel oil brings a realization of 5$/gal, but the flux stock has a cracking value of 8$/gal. This information together with the viscosity and gravity blending numbers and product specifications, appears in the following table. It is desired to operate for maximum profit.
+//
+//--------------------------------------------------------------------
+// Quantity available Value Viscosity Gravity
+// (gal/h) ($/gal) Bl.No. Bl. No
+//Pitch P = 1000-V 0 5 8
+//Visbreaker feed V 0 - -
+//Tar T = 0.8V 0 11 7
+//Flux F = any 8 37 24
+//Fuel oil P + T + F 5 21 min 12 min
+//---------------------------------------------------------------------
+//Suggest a operating condition for maximizing profit.
+
+// The primal model for this example is given as
+//Maximize 5*(P + T + F) - 8*F
+//subject to
+// 5*P + 11*T + 37*F >= 21*(P + T + F)
+// 8*P + 7*T + 24*F >= 12*(P + T + F)
+// P - T/0.8 = 0
+
+// 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;
+//The problem can be execuated using slack variable too.
+//Users are encouraged to explore the use of slack variables.
+
+//Variables = [P T F]; Output = Fuel
+Viscosity = [5 11 37 21];
+Gravity = [8 7 24 12];
+Value = [8 5];
+
+A = [16 10 -16;4 5 -12];
+b = [0 0]';
+Aeq = [0.8 1 0];
+beq = 800;
+lb = [0 0 0];
+ub = [];
+C = -[5 5 -3]';
+Ncons = 3;Nvar = 3;
+//Dual problem
+Adual = -[A;Aeq]';
+bdual = C;
+C = [b;beq];
+// The varable corrosponding to equality constraint is unrestricted
+lb = [0 0 -%inf];
+ub = [];
+[xopt,fopt,exitflag,output,lambda] = linprog(C,Adual,bdual,[],[],lb,ub);
+clc
+
+select exitflag
+ //Display result
+case 0 then
+ mprintf('Optimal Solution Found');
+ input('Press enter to view results');clc;
+ mprintf('Revenue Generated %d\n',fopt);
+ disp(xopt(1:Ncons)','Dual values of the primal constraints');
+
+ for i = 1:Ncons
+ mprintf('\n A cost decrease of %f will occur in the current objective cost value per unit decrease in resource %d ',xopt(i),i)
+ end
+
+case 1 then
+ mprintf('Primal Infeasible')
+case 2 then
+ mprintf('Dual Infeasible')
+case 3
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal')
+case 4
+ mprintf('Solution Abandoned')
+case 5
+ mprintf('Primal objective limit reached')
+else
+ mprintf('Dual objective limit reached')
+end
+
diff --git a/code/linprog/Reddy_Mikks.sce b/code/linprog/Reddy_Mikks.sce
new file mode 100644
index 0000000..805c4e3
--- /dev/null
+++ b/code/linprog/Reddy_Mikks.sce
@@ -0,0 +1,84 @@
+// This is an infeasible problem example. The original problem has been modified by including demand.
+//Ref:H.A. TAHA,"OPERATIONS RESEARCH AN INTRODUCTION",PEARSON-Prentice Hall New Jersey 2007,chapter 2.1
+//Reddy Mikks produces both interior and exterior paints from two raw materials, M1 and M2. The following table provides the basic data of the problem
+//
+// tons of raw material per tons of
+// --------------------------------
+//Raw material Exterior paint Interior paint Maximum daily vailability (tons)
+//-------------------------------------------------------------------------------------------------
+//M1 6 4 24
+//M2 1 2 6
+//-------------------------------------------------------------------------------------------------
+//Profit per ton 5000 4000
+//-------------------------------------------------------------------------------------------------
+//A market survey indicates that the daily demand for interior paint cannot exceed that for exterior paint by more than 1 ton.
+// Also, There should be minimum production of 6 tons combining both paints(this part is changed from the main problem). Reddy Mikks wants to determine the optimum product mix of interior and exterior paints that maximizes the total daily profit.
+
+// 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;
+
+Nproducts = 2;
+Nrawmaterial = 2;
+profit_product =[5000 4000];
+product_Mix = [6 1;4 2];
+Raw_avail = [24 6];
+mprintf('Data received')
+mprintf('Number of product %d',Nproducts);
+mprintf('Number of Raw materials %d',Nrawmaterial);
+FirstR = ['Product 1','Product 2','available'];
+Raw_mat = ['M1','M2'];
+table = [['Raw material' FirstR];[Raw_mat' string(product_Mix') string(Raw_avail') ]];
+disp(table);
+input('Enter to proceed ') ;
+
+for i = 1:Nproducts
+ A(i,:) = product_Mix(:,i)';
+ b(i) = Raw_avail(i);
+end
+
+// market restriction
+A = [A;[-1 1;-1 -1]];
+b = [b(:);1;-6];
+lb = zeros(1,Nproducts);
+ub = [%inf %inf];
+C = profit_product;
+
+[xopt,fopt,exitflag,output,lambda] = linprog(-C',A,b,[],[],lb,ub);
+clc
+select exitflag
+ //Display result
+case 0 then
+ mprintf('Optimal Solution Found')
+ input('Press enter to view results')
+ disp('Paint produced')
+ mprintf('\nExterior paint %d\t interior paint %d',xopt(1),xopt(2));
+case 1 then
+ mprintf('Primal Infeasible')
+ y = A*xopt-b;
+ if sum(y(1:Nrawmaterial))>0
+ mprintf('\nInsufficient raw material');
+ else
+ mprintf('\nCan not meet market constraint');
+ end
+case 2 then
+ mprintf('Dual Infeasible')
+case 3
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal')
+case 4
+ mprintf('Solution Abandoned')
+case 5
+ mprintf('Primal objective limit reached')
+else
+ mprintf('Dual objective limit reached')
+end
+
diff --git a/code/linprog/TOOLCO.sce b/code/linprog/TOOLCO.sce
new file mode 100644
index 0000000..c14ca1b
--- /dev/null
+++ b/code/linprog/TOOLCO.sce
@@ -0,0 +1,82 @@
+//This is an infeasible model
+//Ref:H.A. TAHA,"OPERATIONS RESEARCH AN INTRODUCTION",PEARSON-Prentice Hall New Jersey 2007,chapter 3.5
+//Example:
+//Toolco produces three types of tools, T1,T2 and T3. The tools use two rawmaterials, M1 and M2, according to the data in the following table:
+//--------------------------------------------------------------------
+// Number of units of raw materials per tool
+// -----------------------------------------
+//Raw material T1 T2 T2
+//---------------------------------------------------------------------
+//M1 3 5 6
+//M2 5 3 4
+//---------------------------------------------------------------------
+//The available daily quantites of raw materials M1 and M2 are 1000 units and 1200 units, respectively. The marketing department informed the production manager that according to their research , the daily demand for all three tools must be at least 500 units. will the manufacturing department be able to satisfy the demand?
+
+// Copyright (C) 2015 - 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;
+
+Nraw = 2;
+Nproduct = 3;
+Product_mix = [3 5 6;5 3 4];
+Raw_avail = [1000 1200];
+Demand = [500 500 500];
+
+mprintf('Data received')
+mprintf('Number of product %d',Nproduct);
+mprintf('Number of Raw materials %d',Nraw);
+FirstR = ['Product 1','Product 2','Product 3','Available']
+Raw_mat = ['M1','M2'];
+
+table = [['Raw materisl' FirstR];[Raw_mat' string(Product_mix) string(Raw_avail') ]];
+disp(table);
+input('Enter to proceed ');
+
+for i = 1:Nraw
+ A(i,:) = Product_mix(i,:);
+ b(i) = Raw_avail(i);
+end
+A = [A;-eye(Nproduct,Nproduct)];
+b = [b;-500*ones(Nproduct,1)];
+
+lb = zeros(1,Nproduct);//Redudant due to demand
+ub = [];
+C = ones(1,Nproduct);
+[xopt,fopt,exitflag,output,lambda] = linprog(-C',A,b,[],[],lb,ub);
+
+clc
+select exitflag
+ //Display result
+case 0 then
+ mprintf('Optimal Solution Found')
+ input('Press enter to view results')
+ disp('Paint produced')
+ mprintf('\nExterior paint %d\t interior paint %d',xopt(1),xopt(2));
+case 1 then
+ mprintf('Primal Infeasible\n')
+ y = A*xopt-b;
+ if sum(y(1:Nraw))>0
+ mprintf('\nInsufficient raw material');
+ else
+ mprintf('\nCan not meet market demand');
+ end
+case 2 then
+ mprintf('Dual Infeasible')
+case 3
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal')
+case 4
+ mprintf('Solution Abandoned')
+case 5
+ mprintf('Primal objective limit reached')
+else
+ mprintf('Dual objective limit reached')
+end
+
diff --git a/code/linprog/TOYCO Dual.sce b/code/linprog/TOYCO Dual.sce
new file mode 100644
index 0000000..8fdba97
--- /dev/null
+++ b/code/linprog/TOYCO Dual.sce
@@ -0,0 +1,88 @@
+//This example shows the use of duality theorm to determine the reduced cost associated with the non-negative constraints
+//Ref:H.A. TAHA,"OPERATIONS RESEARCH AN INTRODUCTION",PEARSON-Prentice Hall New Jersey 2007,chapter 2.1
+//Example:
+//TOYCO assembels three types of toys-trains,trucks and cars using three operations. the daily limits on the available times for the three operations are 430,460 and 420 minutes respectively. The revenues per unit of toy train, truck and car are $3.$2 and $5 respectively. The corrosponding times per train and per car are (2,0,4) and (1,2,0) minutes ( a zero time indicates that the operation is not used). determine the . Determine the shadow pricing and how optima changes with resource change.
+
+//The primal model of the above problem is
+//maximize z = 3*x1 + 2*x2 + 5*x3
+//Subject to
+//x1 + 2*x2 + x3 <=430
+//3*x1 +2*x3 <=460
+//x1 + 4*x2 <=420
+//x1,x2,x3 >= 0
+
+// 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;
+
+Nproducts = 3;
+Noperations = 3;
+revenue = [3 2 5];
+resource = [430 460 420];
+Assembly_time = [1 3 1;2 0 4;1 2 0];
+mprintf('Data Received')
+prodNam = ['Trains','Trucks','cars']';
+Operation = ['Operation 1','Operation 2','Operation 3'];
+table = [['Products',Operation,'Revenue'];[prodNam,string(Assembly_time),string(revenue)'];['Resource',string(resource),'']];
+disp(table)
+input('Press enter to proceed ')
+//primal model constraints
+A = Assembly_time';
+b = resource';
+
+//dual constraints
+dualA = -A';
+[Ncons,Nvars] = size(A);
+for i = 1:Nvars
+ A1(i,i) = 1;
+end
+dualAeq = [dualA A1];
+dualbeq = -revenue';
+C = [b;zeros(Nvars,1)];
+lb= zeros(1,2*Nvars);
+ub = [];
+[xopt,fopt,exitflag,output,lambda] = linprog(C,[],[],dualAeq,dualbeq,lb,[]);
+
+clc;
+select exitflag
+ //Display result
+case 0 then
+ mprintf('Optimal Solution Found');
+ input('Press enter to view results');clc;
+ mprintf('Revenue Generated %d\n',fopt);
+ disp(xopt(1:Ncons)','Dual values of the primal constraints');
+ disp(xopt(Ncons+1:Ncons+Nvars)','Dual values of each primal variable');
+ for i = 1:Ncons
+ mprintf('\n A cost decrease of %d will occur in the current objective cost value per unit decrease in resource %d ',xopt(i),i)
+ end
+ for i = 1:Nvars
+ mprintf('\n A cost decrease of %d will occur in the current objective cost value per unit increase in variable x%d ',xopt(Ncons+i),i)
+ end
+
+case 1 then
+ mprintf('Primal Infeasible')
+case 2 then
+ mprintf('Dual Infeasible')
+case 3
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal')
+case 4
+ mprintf('Solution Abandoned')
+case 5
+ mprintf('Primal objective limit reached')
+else
+ mprintf('Dual objective limit reached')
+end
+
+
+
+
+
diff --git a/code/linprog/TOYCO model.sce b/code/linprog/TOYCO model.sce
new file mode 100644
index 0000000..b6ea553
--- /dev/null
+++ b/code/linprog/TOYCO model.sce
@@ -0,0 +1,67 @@
+//This example shows the use of slack or surplus variables and their sigificance on the model
+//Ref:H.A. TAHA,"OPERATIONS RESEARCH AN INTRODUCTION",PEARSON-Prentice Hall New Jersey 2007,chapter 2.1
+//Example:
+//TOYCO assembels three types of toys-trains,trucks and cars using three operations. the daily limits on the available times for the three operations are 430,460 and 420 minutes respectively. The revenues per unit of toy train, truck and car are $3.$2 and $5 respectively. The corrosponding times per train and per car are (2,0,4) and (1,2,0) minutes ( a zero time indicates that the operation is not used). determine the optimum production rate to maximize profit. Check any surplus resource availability.
+
+// 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;
+
+Nproducts = 3;
+Noperations = 3;
+revenue = [3 2 5];
+resource = [430 460 420];
+Assembly_time = [1 3 1;2 0 4;1 2 0];
+mprintf('Data Received')
+prodNam = ['Trains','Trucks','cars']';
+Operation = ['Operation 1','Operation 2','Operation 3'];
+table = [['Products',Operation,'Revenue'];[prodNam,string(Assembly_time),string(revenue)'];['Resource',string(resource),'']];
+disp(table)
+input('Press enter to proceed ')
+
+NequalCon = Noperations;
+Aeq = zeros(NequalCon,Nproducts+Noperations);
+for i = 1:Noperations
+ Aeq(i,1:Nproducts) = Assembly_time(:,i)';
+ beq(i) = resource(i);
+ Aeq(i,i+Nproducts) = 1;
+end
+
+C = -[revenue zeros(1,Nproducts)];
+lb= zeros(1,2*Noperations);
+ub = [];
+
+[xopt,fopt,exitflag,output,lambda] = linprog(C',[],[],Aeq,beq,lb,[]);
+clc;
+select exitflag
+ //Display result
+case 0 then
+ mprintf('Optimal Solution Found');
+ input('Press enter to view results');clc;
+ mprintf('Revenue Generated %d\n',-fopt);
+ disp(xopt(1:Nproducts)','Optimal Number of Trains Trucks and Cars');
+ disp(xopt(Nproducts+1:Nproducts+NequalCon)','Surplus resource available among 3 resources');
+ beq = beq-xopt(Nproducts+1:Nproducts+NequalCon);
+ disp(beq','Minimum Constraint boundry for Current optima');
+case 1 then
+ mprintf('Primal Infeasible')
+case 2 then
+ mprintf('Dual Infeasible')
+case 3
+ mprintf('Maximum Number of Iterations Exceeded. Output may not be optimal')
+case 4
+ mprintf('Solution Abandoned')
+case 5
+ mprintf('Primal objective limit reached')
+else
+ mprintf('Dual objective limit reached')
+end
generated by cgit (git 2.34.1) at 2024-12-19 10:18:36 +0000