From 69460c03b8b53068d60fd08d3180efc91e627603 Mon Sep 17 00:00:00 2001 From: RemyaDebasis Date: Mon, 23 Jul 2018 20:01:22 +0530 Subject: added code files --- code/Symphonymat/Factory planing.sce | 143 ++++++++++++++++++++++++++++++ code/Symphonymat/Mining.sce | 154 +++++++++++++++++++++++++++++++++ code/Symphonymat/factory.xls | Bin 0 -> 29184 bytes code/fminbnd/Chichinadze.sce | 64 ++++++++++++++ code/fminbnd/Cola.sce | 89 +++++++++++++++++++ code/fmincon/Tankprob.sci | 95 ++++++++++++++++++++ code/fmincon/spring.sci | 120 +++++++++++++++++++++++++ code/fminsearch/Brownsfunc.sci | 73 ++++++++++++++++ code/fminsearch/FletcherPowell.sce | 65 ++++++++++++++ code/fminunc/Freudenstein and Roth.sci | 75 ++++++++++++++++ code/fsolve/Redlich-Kwong.sce | 69 +++++++++++++++ code/fsolve/SLEs.sce | 50 +++++++++++ code/intfmincon/ProcessSel.sce | 148 +++++++++++++++++++++++++++++++ code/intqpipopt/SpringCart.sce | 68 +++++++++++++++ code/linprog/Blending_problem.sce | 135 +++++++++++++++++++++++++++++ code/linprog/Employee Scheduling.sce | 83 ++++++++++++++++++ code/linprog/FuelOil.sce | 86 ++++++++++++++++++ code/linprog/Reddy_Mikks.sce | 84 ++++++++++++++++++ code/linprog/TOOLCO.sce | 82 ++++++++++++++++++ code/linprog/TOYCO Dual.sce | 88 +++++++++++++++++++ code/linprog/TOYCO model.sce | 67 ++++++++++++++ 21 files changed, 1838 insertions(+) create mode 100644 code/Symphonymat/Factory planing.sce create mode 100644 code/Symphonymat/Mining.sce create mode 100644 code/Symphonymat/factory.xls create mode 100644 code/fminbnd/Chichinadze.sce create mode 100644 code/fminbnd/Cola.sce create mode 100644 code/fmincon/Tankprob.sci create mode 100644 code/fmincon/spring.sci create mode 100644 code/fminsearch/Brownsfunc.sci create mode 100644 code/fminsearch/FletcherPowell.sce create mode 100644 code/fminunc/Freudenstein and Roth.sci create mode 100644 code/fsolve/Redlich-Kwong.sce create mode 100644 code/fsolve/SLEs.sce create mode 100644 code/intfmincon/ProcessSel.sce create mode 100644 code/intqpipopt/SpringCart.sce create mode 100644 code/linprog/Blending_problem.sce create mode 100644 code/linprog/Employee Scheduling.sce create mode 100644 code/linprog/FuelOil.sce create mode 100644 code/linprog/Reddy_Mikks.sce create mode 100644 code/linprog/TOOLCO.sce create mode 100644 code/linprog/TOYCO Dual.sce create mode 100644 code/linprog/TOYCO model.sce 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 Binary files /dev/null and b/code/Symphonymat/factory.xls 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 j2 + 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= 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 -- cgit