diff options
Diffstat (limited to 'code')
24 files changed, 1806 insertions, 0 deletions
diff --git a/code/Symphonymat/Assignment_Problem.sce b/code/Symphonymat/Assignment_Problem.sce new file mode 100644 index 0000000..419d9ee --- /dev/null +++ b/code/Symphonymat/Assignment_Problem.sce @@ -0,0 +1,79 @@ +//Reference :H.A. Eiselt and C.-L.SAndblom,"Linear Programming and its Applications",Springer-Verlag Berlin Heidelberg 2007,chapter 2.9 + +//The public works department of a region has issued contracts for five new projects. Each of the three contractors who have shown interest in the projects has submitted bids for those projects. The bids given by contactors are as given below +// ===================================================================== +// Project1 Project2 Project3 Project4 Project5 +// --------------------------------------------------------------------- +// Contractor1 20 - 10 9 15 +// Contractor2 18 12 13 8 16 +// Contractor3 - 11 12 7 17 +// ===================================================================== +//The resource requirement of the contractors for each projects and the total resources available with them are as given +// ============================================================================= +// Project1 Project2 Project3 Project4 Project5 Available Resources +// ----------------------------------------------------------------------------- +// Contractor1 70 - 40 30 50 120 +// Contractor2 65 45 40 35 55 100 +// Contractor3 - 45 35 40 50 70 +// ============================================================================ +//The objective of the contractors is to minimize the total cost +clc; +nProjects = 5; +nContractors = 3; +bids = [20 10000 10 9 15;18 12 13 8 16; 10000 11 12 7 17]; +reqResource = [70 10000 40 30 50;65 40 40 35 55; 10000 45 35 40 50]; +availResource = [120; 100; 70]; + +nVar = nProjects*nContractors; // Calculated the dimension of the problem +IntCon = 1:nVar; // Indicating the integer variables +lb = zeros(1,nVar); +ub = ones(1,nVar); + +// Linear constraints +// Linear equality constraints + +beq = ones(nProjects,1); +for i = 1:nProjects + Aeq(i,i:nProjects:nVar) = 1; +end + +// Linear inequality constraints +b = availResource; +for j = 1:nContractors + index = (j-1)*nProjects+1:j*nProjects; + A(j,index) = reqResource(j,:); +end + +// Objective function +for j = 1:nContractors + index = (j-1)*nProjects+1:j*nProjects; + Cost(index,1) = bids(j,:)'; +end + +options = list("time_limit", 2500); +[xopt,fopt,status,output] = symphonymat(Cost,IntCon,A,b,Aeq,beq,lb,ub,options); + +// Result representation +select status +case 227 + disp(" Optimal Solution Found") +case 228 + disp("Maximum CPU Time exceeded") +case 229 + disp("Maximum Number of Node Limit Exceeded") +case 230 + disp("Maximum Number of Iterations Limit Exceeded.") +end + +for j = 1:nContractors + index = (j-1)*nProjects+1:j*nProjects; + Contractor(j).Project = string(find(xopt(index)==1)); +end +ResolurceUtilized = ((A*xopt)./b)*100; + +for j = 1:nContractors + disp(strcat(["Projects assigned to contractor",string(j)," : ", Contractor(j).Project],' ')); + disp(strcat(["Resource utilized by contractor ",string(j)," : ",string(ResolurceUtilized(j)),"%"])); +end +disp(strcat(["Total cost : ",string(fopt)])) + diff --git a/code/fgoalattain/WeldedBealGoalProgramming.sce b/code/fgoalattain/WeldedBealGoalProgramming.sce new file mode 100644 index 0000000..6bf1242 --- /dev/null +++ b/code/fgoalattain/WeldedBealGoalProgramming.sce @@ -0,0 +1,89 @@ +//Reference: K. Deb,"Solving Goal Programming Problems Using Multi-Objective Genetic Algorithms",Proceeedings of the 1999 Congress on Evolutionary C omputation CEC-99,USA, 1999, p. 77-84 +// goal f1(x) = (1.10471*h^2*l + 0.04811*t*b*(14+l))<=5; +// goal f2(x) = (2.1952/((t^3)*b))<=0.001; +// subjected to +// g1(x) = 13,600 - tau >= 0 +// g2(x) = 30,000 - sigma >= 0 +// g3(x) = b - h >= 0 +// g4(x) = Pc(x) - 6000 >= 0 +// 0.125 <= h,b <= 5 and 0.1 <= l,t <= 10 +//===================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== +clc; +// Objective functions +function f = ObjectiveFunction(X) + f(1) = (1.10471*X(1)^2*X(2) + 0.04811*X(3)*X(4)*(14+X(2))); + f(2) = (2.1952/((X(3)^3)*X(4))); +endfunction +// Non linear constraints +function [C,Ceq] = NLconstraints(X) + P = 6000; L = 14; E = 3*10^7; G = 12*10^6; tauMax = 13600; sigmaMax = 30000; + + M = P*(L+(X(2)/2)); + R = sqrt((X(2)^2/4) + ((X(1) + X(3))/2)^2); + J = 2*0.7071068*X(1)*X(2)*((X(2)^2/12)+((X(1)+X(3))/2)^2); + sigma = (6*P*L)/(X(3)^2*X(4)); + Pc1 = (4.013*(sqrt(E*G*(X(3)^2*X(4)^6)/36)))/(L^2); + Pc2 = 1-(X(3)/(2*L))*sqrt(E/(4*G)); + Pc = Pc1*Pc2; + tauPrime = P/(sqrt(2)*X(1)*X(2)); + tauDprime = (M*R)/J; + tau = sqrt(tauPrime^2 + 2*tauPrime*tauDprime*(X(2)/(2*R))+tauDprime^2); + + C(1) = tau - tauMax; + C(2) = sigma - sigmaMax; + C(3) = P - Pc; + C = C'; + Ceq = []; +endfunction +// Linear inequality constraints +A = [1 0 0 -1]; +b = 0; + +nObj = 2; + +Aeq = []; beq = []; +lb = [0.125 0.1 0.1 0.125]; +ub = [5 10 10 5]; +nVar = length(lb); + +// Initial guess to the solver +x0 = lb + rand(1,nVar).*(ub-lb); + +//Specifying the goal and the weights +goal=[5 0.001]; +weight = [0.75 0.25]; + +options = list("MaxIter",10000,"CpuTime",10000) +// Calling fgoalattain +[x,fval,attainfactor,exitflag,output,lambda]=fgoalattain(ObjectiveFunction,x0,goal,weight,A,b,Aeq,beq,lb,ub,NLconstraints,options) + +// Result representation +clc; +select exitflag +case 0 + disp("Optimal Solution Found") + disp(x0,"Initial guess given to the solver") + disp(x',"The optimum solution obtained") + disp(fval',"The optimum value of the objective functions") +case 1 + disp(" Maximum Number of Iterations Exceeded. Output may not be optimal") +case 2 + disp("Maximum amount of CPU Time exceeded. Output may not be optimal") +case 3 + disp("Stop at Tiny Step") +case 4 + disp("Solved To Acceptable Level") +case 5 + disp("Converged to a point of local infeasibility") +end + diff --git a/code/fminbnd/LevyFunction.sce b/code/fminbnd/LevyFunction.sce new file mode 100644 index 0000000..2a0d4ee --- /dev/null +++ b/code/fminbnd/LevyFunction.sce @@ -0,0 +1,69 @@ +//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 + +// n-dimensional Levy Function +// Dimension,D = [4 5 10 50 100]; +// Domain = |Xi| <= 10 where i = 1,2,...,D +// Global minimum is zero at Xi = 1. +//====================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//====================================================================== +clc; +// Objective function +function f = LevyFunction(X) + d = length(X); + w1 = 1 + (X(1)-1)/4; + wi = 1 + (X(1:d-1)-1)./4; + wd = 1 + (X(d)-1)/4; + Term1 = (sin(%pi.*w1)).^2; + Term2 = sum(((wi-1).^2).*(1+10*(sin(%pi.*wi + 1)).^2)); + Term3 = ((wd-1).^2)*(1+(sin(2*%pi*X(d))).^2); + f = Term1 + Term2 + Term3; +endfunction +// Dimension of the problem +D = 5; +// Bounds of the problem +LB = -10*ones(1,D); +UB = 10*ones(1,D); +// Calling the solver +options=list("TolX",0) +[xopt,fopt,exitflag,output,lambda] = fminbnd(LevyFunction,LB,UB,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") + disp(fopt,"The objective function value is") +case 2 + disp("Maximum CPU Time exceeded. Output may not be optimal") + disp(xopt,"The solution obtained") + disp(fopt,"The objective function value is") +case 3 + disp("Stop at Tiny Step") + disp(xopt,"The solution obtained") + disp(fopt,"The objective function value is") +case 4 + disp("Solved To Acceptable Level") + disp(xopt,"The solution obtained") + disp(fopt,"The objective function value is") +case 5 + disp("Converged to a point of local infeasibility") + disp(xopt,"The solution obtained") + disp(fopt,"The objective function value is") +end +disp(output) diff --git a/code/fminbnd/MichalewiczFunction.sce b/code/fminbnd/MichalewiczFunction.sce new file mode 100644 index 0000000..b700973 --- /dev/null +++ b/code/fminbnd/MichalewiczFunction.sce @@ -0,0 +1,29 @@ +//Reference: Michalewicz, Z.: Genetic Algorithms + Data Structures = Evolution Programs. Berlin, Heidelberg, New York: Springer-Verlag, 1992 + +// MICHALEWICZ FUNCTION +clc; +function f = ObjectiveFunction(X) + m = 10; + nVar = length(X); + d = length(X); + f = 0; + for n = 1:nVar + f = f - sin(X(n))*((sin((n*X(n)^2)/%pi))^(2*m)); + end + + f = -f; +endfunction + +nVar = 2; +lb = zeros(1,nVar); +ub = %pi*ones(1,nVar); +[xopt,fopt,exitflag,output,lambda] = fminbnd(ObjectiveFunction,lb,ub) + +disp (nVar) +disp(fopt) +disp(xopt') + + + + + diff --git a/code/fmincon/Continuous through circulation dryer.sce b/code/fmincon/Continuous through circulation dryer.sce new file mode 100644 index 0000000..bbf38ab --- /dev/null +++ b/code/fmincon/Continuous through circulation dryer.sce @@ -0,0 +1,66 @@ +//Reference: Rein Luus and Taina . I. Jaakola, Optimization of Nonlinear Functions Subject to Equality Constraints. Judicious Use of Elementary Calculus and Random Numbers,Ind. Eng. Chem. Process Des. Develop., Vol. 12, No. 3, 1973 + +// Optimization of a continuous through circulation dyer +// Maximize the production rate by considering the fluid velocity (X1) and bed length (X2) as the design parameters subjected to power constraint and the moisture content distribution constraint + +//====================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//====================================================================== +clc; +// Objective function +function f = ObjectiveFunction(X) + f = -(0.0064*X(1)*(1 - exp(-0.184*(X(1)^0.3)*X(2)))); +endfunction +// Non linear equality constraints +function [C,Ceq] = NLConstraints(X) + Ceq(1,1) = (3000 + X(1))*(X(1)^2)*X(2) - 1.2D13; + Ceq(1,2) = exp(0.184*X(1)^0.3*X(2)) - ln(4.1); + C = []; +endfunction +// Bounds of the problem +lb = [0 0]; +ub = []; +// Initial guess +x0 = [1000 20]; +// Calling the solver +[xopt,fopt,exitflag,output,lambda] = fmincon(ObjectiveFunction,x0,[],[],[],[],lb,[],NLConstraints) + +// Result representation +select exitflag +case 0 + disp("Optimal Solution Found") + disp(x0,"Initial guess given to the solver") + disp(xopt',"The optimum solution obtained") + disp(fopt,"The optimum value of the objective function") +case 1 + disp(" Maximum Number of Iterations Exceeded. Output may not be optimal") + disp(x0,"Initial guess given to the solver") + disp(xopt',"The solution obtained") + disp(fopt,"The value of the objective function") +case 2 + disp("Maximum amount of CPU Time exceeded. Output may not be optimal") + disp(x0,"Initial guess given to the solver") + disp(xopt',"The solution obtained") + disp(fopt,"The value of the objective function") +case 3 + disp("Stop at Tiny Step") + disp(x0,"Initial guess given to the solver") + disp(xopt',"The solution obtained") + disp(fopt,"The value of the objective function") +case 4 + disp("Solved To Acceptable Level") + disp(x0,"Initial guess given to the solver") + disp(xopt',"The solution obtained") + disp(fopt,"The value of the objective function") +case 5 + disp("Converged to a point of local infeasibility") +end + diff --git a/code/fmincon/Welded Beam Design.sce b/code/fmincon/Welded Beam Design.sce new file mode 100644 index 0000000..b520f31 --- /dev/null +++ b/code/fmincon/Welded Beam Design.sce @@ -0,0 +1,86 @@ +//Reference: K. M. Ragsdell and D. T. Phillips,"Optimal Design of a Class of Welded Structures Using Geometric Programming",ASME Journal of Engineering for Industry,Vol.98, pp 1021-1025, 1976 +//A welded beam is designed for minimum cost subject to constraints on shear stress,bending stress in the beam,buckling load on the bar,end deflectionof the beam and the bound constraints. There are four design variables such as Weld thickness(inch), weld length (inch), bar thickness (inch), bar breadth(inch). +//===================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== +clc; +// Objective function +function f = ObjectiveFunction(X) + f = 1.10471*X(1)^2*X(2) + 0.04811*X(3)*X(4)*(14+X(2)); +endfunction +// Non linear equality and inequality constraints +function [C,Ceq] = NLconstraints(X) + P = 6000; L = 14; E = 3*10^7; G = 12*10^6; tauMax = 13600; + sigmaMax = 30000; deltaMax = 0.25; + + M = P*(L+(X(2)/2)); + R = sqrt((X(2)^2/4) + ((X(1) + X(3))/2)^2); + J = 2*0.7071068*X(1)*X(2)*((X(2)^2/12)+((X(1)+X(3))/2)^2); + sigma = (6*P*L)/(X(3)^2*X(4)); + delta = (4*P*L^3)/(E*X(3)^3*X(4)); + Pc1 = (4.013*(sqrt(E*G*(X(3)^2*X(4)^6)/36)))/(L^2); + Pc2 = 1-(X(3)/(2*L))*sqrt(E/(4*G)); + Pc = Pc1*Pc2; + tauPrime = P/(sqrt(2)*X(1)*X(2)); + tauDprime = (M*R)/J; + tau = sqrt(tauPrime^2 + 2*tauPrime*tauDprime*(X(2)/(2*R))+tauDprime^2); + + C(1) = tau - tauMax; + C(2) = sigma - sigmaMax; + C(3) = 0.10471*X(1)^2 + 0.04811*X(3)*X(4)*(14+X(2)) - 5; + C(4) = delta - deltaMax; + C(5) = P - Pc; + C = C'; + Ceq = []; +endfunction +// Linear inequality constraints +A = [1 0 0 -1;-1 0 0 0]; +b = [0 -0.125]'; +// Bounds of the problem +lb = [0.1 0.1 0.1 0.1]; +ub = [2 10 10 2]; +// Initial guess +x0 = rand(1,4).*(ub-lb); +// Design parameters of the problem +designParameters = {'Weld thickness(inch)','weld length (inch)','bar thickness (inch)','bar breadth(inch)'} +inGuess = [designParameters; string(x0)] +disp(inGuess,"Initial guess given to the solver") +input("Press enter to proceed: ") +// Calling the solver function +[xopt,fopt,exitflag,output,lambda] = fmincon(ObjectiveFunction,x0,A,b,[],[],lb,ub,NLconstraints) + +// Result representation +clc; +optSol = [designParameters; string(xopt')] +select exitflag +case 0 + disp(optSol,"The optimum solution obtained") + disp(fopt,"The minimum cost for the weldment assembly is") +case 1 + disp(" Maximum Number of Iterations Exceeded. Output may not be optimal") + disp(optSol,"The solution obtained") + disp(fopt,"The minimum cost for the weldment assembly is") +case 2 + disp("Maximum amount of CPU Time exceeded. Output may not be optimal") + disp(optSol,"The solution obtained") + disp(fopt,"The minimum cost for the weldment assembly is") +case 3 + disp("Stop at Tiny Step") + disp(optSol,"The solution obtained") + disp(fopt,"The minimum cost for the weldment assembly is") +case 4 + disp("Solved To Acceptable Level") + disp(optSol,"The solution obtained") + disp(fopt,"The minimum cost for the weldment assembly is") +case 5 + disp("Converged to a point of local infeasibility") +end +disp(output) diff --git a/code/fminimax/Davidson2Problem.sce b/code/fminimax/Davidson2Problem.sce new file mode 100644 index 0000000..e00bdc2 --- /dev/null +++ b/code/fminimax/Davidson2Problem.sce @@ -0,0 +1,52 @@ +// J.Hald, K Madsen, Reference:Combined LP and quasi-Newton methos for minimax optimization, Mathematical Programming, Springer, 1981 +// Min F = max|fi(x)| +// fi(X) = (X1 + X2*ti - exp(ti))^2 + (X3 + X4*sin(ti) - cos(ti))^2 where i varies from 1 to 20 +// ti = 0.2i +// Global optima: X* = [-12.24368 14.02180 -0.45151 -0.01052]; F* = 115.70644; +//====================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//====================================================================== +clc; +// Objective function +function f = ObjectiveFunction(X) + m = 20; + for i = 1:m + t(i) = 0.2*i; + f(i) = abs((X(1) + X(2)*t(i) - exp(t(i)))^2 + (X(3) + X(4)*sin(t(i)) - cos(t(i)))^2); + end +endfunction +// Initial guess to the solver +x0 = [-10 5 1 -2] + +// Run fminimax +options= list("MaxIter", 10000); +[x,fval,maxfval,exitflag,output,lambda] = fminimax(ObjectiveFunction,x0,[],[],[],[],[],[],[],options); + +// Result representation +clc; +select exitflag +case 0 + disp("Optimal Solution Found") + disp(x0,"The initial guess given to the solver") + disp(x',"The optimal solution determined by the solver") + disp(maxfval,"The optimum value of the objective function") +case 1 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 2 + disp("Maximum amount of CPU Time exceeded. Output may not be optimal") +case 3 + disp("Stop at Tiny Step") +case 4 + disp("Solved To Acceptable Level") +case 5 + disp("Converged to a point of local infeasibility") +end + diff --git a/code/fminimax/NL_MiniMax_Problem.sce b/code/fminimax/NL_MiniMax_Problem.sce new file mode 100644 index 0000000..759f56b --- /dev/null +++ b/code/fminimax/NL_MiniMax_Problem.sce @@ -0,0 +1,63 @@ +//C. Charalambous and A.R. Conn, "An efficient method to solve the minimax problem directly",SIAM Journal on Numerical Analysis 15 (1978) 162-187 + +//====================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//====================================================================== +clc +// Objective function +function f = ObjectiveFunction(X) + g(1) = 0; + g(2) = -2*X(1)^2 - 3*X(2)^4 - X(3) - 4*X(4)^2 - 5*X(5) + 127; + g(3) = -7*X(1) - 3*X(2) - 10*X(3)^2 - X(4) + X(5) + 282; + g(4) = -23*X(1) - X(2)^2 - 6*X(6)^2 + 8*X(7) + 196; + g(5) = -4*X(1)^2 - X(2)^2 + 3*X(1)*X(2) - 2*X(3)^2 - 5*X(6) + 11*X(7); + f(1) = (X(1) - 10)^2 + 5*(X(2) - 12)^2 + X(3)^4 + 3*(X(4) - 11)^2 + 10*X(5)^6 + 7*X(6)^2 + X(7)^4 - 4*X(6)*X(7) - 10*X(6) - 8*X(7); + for i = 2:5 + f(i) = f(1) - 10*g(i); + end + f1(2:5) = f(1) + 10*g(2:5); +endfunction +// Non linear inequality constraints +function [C, Ceq] = NLConstraints(X) + C(1) = -(-2*X(1)^2 - 3*X(2)^4 - X(3) - 4*X(4)^2 - 5*X(5) + 127); + C(2) = -(-7*X(1) - 3*X(2) - 10*X(3)^2 - X(4) + X(5) + 282); + C(3) = -(-23*X(1) - X(2)^2 - 6*X(6)^2 + 8*X(7) + 196); + C(4) = -(-4*X(1)^2 - X(2)^2 + 3*X(1)*X(2) - 2*X(3)^2 - 5*X(6) + 11*X(7)); + C = C'; + Ceq = []; +endfunction + +A =[]; b = []; Aeq = []; beq = [];lb = [];ub = []; +// Initial guess given to the solver +x0 = [2.33050, 1.95137, -0.47754, 4.36573, -0.62449, 1.03813, 1.59423]; +// Solving the solvers +options = list("maxiter", 4000, "cputime", [6000]); +[x,fval,maxfval,exitflag,output,lambda] = fminimax(ObjectiveFunction,x0,A,b,Aeq,beq,lb,ub,NLConstraints,options); + +// Result representation +clc +select exitflag +case 0 + disp("Optimal Solution Found") + disp(x0,"The initial guess given to the solver") + disp(x',"The optimal solution determined by the solver") + disp(maxfval,"The optimum value of the objective function") +case 1 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 2 + disp("Maximum amount of CPU Time exceeded. Output may not be optimal") +case 3 + disp("Stop at Tiny Step") +case 4 + disp("Solved To Acceptable Level") +case 5 + disp("Converged to a point of local infeasibility") +end diff --git a/code/fminsearch/MinOilCost.sce b/code/fminsearch/MinOilCost.sce new file mode 100644 index 0000000..8ba1c94 --- /dev/null +++ b/code/fminsearch/MinOilCost.sce @@ -0,0 +1,81 @@ +//Reference: Edgar, Himmelblau and Lasdon,"Optimization of Chemical Processes",2nd Ed,McGraw-Hill Chemical Engineering Series,chapter 6 + +//The cost of refined oil when shipped via the Malacca Straits to Japan in dollars per kiloliter was given as the linear sum of the crude oil cost, the insurance, customs, freight cost for the oil, loading and unloading cost, sea berth cost, submarine pipe cost, storage cost, tank areacost, refining cost, and freight cost of products. Compute the minimum cost of oil and the optimum tanker size t and refinery size. + +//====================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//====================================================================== +clc; + +function f = ObjectiveFunction(X) + t = X(1); q = X(2); + Cc = 12.5; // Crude oil price ($/kL) + Ci = 0.5; // insurance cost ($/kL) + Cx = 0.9; // customs cost ($/kL) + a = 0.2; // anual fixed charges,fraction + p = 7000; // land prices ($/square meter) + n = 2; // number of ports + i = 0.1; // interest rate + Term1 = (52.47*q*360); + Term2 = (n*t+1.2*q); + f = Cc + Ci + Cx + (2.09 *10^4*t^-0.3017)/360 +... + (1.064*10^6*a*t^0.4925)/Term1+... + (4.242*10^4*a*t^0.7952 + 1.813*i*p*Term2^0.861)/Term1... + +(4.25*10^3*a*Term2)/Term1 + (5.042*10^3*q^-0.1899)/360 +... + (0.1049*q^0.671)/360; +endfunction + + +function stop=outfun(x, optimValues, state) + subplot(1,2,1) + plot(optimValues.iteration,optimValues.fval,'rp'); + xlabel('Iteration');ylabel('fval') + + subplot(1,2,2) + plot(optimValues.iteration,x(1),'b*'); + plot(optimValues.iteration,x(2),'g*'); + legend(['Tanker size','Refinery size']) + set(gca(),"auto_clear","off") + xlabel('Iteration');ylabel('variable values') + + stop = %f +endfunction +designParameters = {"Tanker size","Refinery size"} +X0 = [15000 20000]; +designParameter = {'Tanker size(kL)','Refinery capacity(bbl/day)'}; +intGuess = [designParameter;string(X0)]; +disp(intGuess,"Initial guess given to the solver") +input('Press enter to proceed: ') +Parameter.X_Tol = 0; +Parameter.F_Tol = 0; +Parameter.maxFE = 1000; +Parameter.maxIt = 100; +opt = optimset ("TolX",Parameter.X_Tol,"TolFun",Parameter.F_Tol,"MaxFunEvals",Parameter.maxFE,"MaxIter",Parameter.maxIt,"OutputFcn",outfun); + +[x,fval,exitflag,output] = fminsearch(ObjectiveFunction,X0,opt) + +clc; +select exitflag +case -1 + mprintf('The maximum number of iterations has been reached \n') + mprintf('Function Count: %d ',output.funcCount) +case 0 + mprintf('The maximum number of function evaluations has been reached \n') + mprintf(' Iteration Count: %d ',output.iterations) + +case 1 + mprintf('The tolerance on the simplex size and function value delta has been reached \n') + mprintf('Function Count: %d ',output.funcCount) + mprintf('Iteration Count: %d ',output.iterations) +end +optSol = [designParameter;string(x)]; +disp(optSol,"The optimum solution obtained") +disp(fval,"The minimum cost of the oil is") diff --git a/code/fminsearch/PowellsBadlyScaledFunction.sci b/code/fminsearch/PowellsBadlyScaledFunction.sci new file mode 100644 index 0000000..531ac2d --- /dev/null +++ b/code/fminsearch/PowellsBadlyScaledFunction.sci @@ -0,0 +1,48 @@ +// This is an example for unconstrained multivariable problem +// Name of the function : Powell's Badly Scaled function +// f(x1,x2) = (10000 x1.x2 - 1)^2 + [exp(-x1) + exp(-x2) - 1.0001]^2 + +// Reference: M.J.D Powell, A hybrid method for non-linear equations, pp.87-114 in numerical methods for non-linear algebraic equations,P.Rabinowitz,Ed.,Gorden and Breach,Newyork,1970 +//===================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== +clc; + +function f = PowellsBadlyScaled(X) + f = (10000*X(1)*X(2) -1)^2 + (exp(-X(1)) + exp(-X(2)) - 1.0001)^2; +endfunction + +X0 = [0 0]; +Parameter.X_Tol = 1.e-16; +Parameter.F_Tol = 1.e-16; +Parameter.maxFE = 1000; +Parameter.maxIt = 400; +mprintf("The values set for the configurable options are as below") +disp(Parameter); +input("Press enter to proceed ") + +opt = optimset ("TolX",Parameter.X_Tol,"TolFun",Parameter.F_Tol,"MaxFunEvals",Parameter.maxFE,"MaxIter",Parameter.maxIt,"PlotFcns" , optimplotfval); + +[x,fval,exitflag,output] = fminsearch(PowellsBadlyScaled,X0,opt) + +clc +select exitflag +case -1 + mprintf('The maximum number of iterations has been reached \n') +case 0 + mprintf('The maximum number of function evaluations has been reached \n') +case 1 + mprintf('The tolerance on the simplex size and function value delta has been reached \n') +end + +disp(x,"The optimal solution is") +disp(fval,"The optimum value of the objective function is") +disp(output) diff --git a/code/fminunc/SStemperature.sce b/code/fminunc/SStemperature.sce new file mode 100644 index 0000000..cdbd08d --- /dev/null +++ b/code/fminunc/SStemperature.sce @@ -0,0 +1,53 @@ +//Reference: S.S. Rao,Engineering Optimization Theory and Practice, 3rd enlarged edition, New age international publishers,2011,chapter 6 +// The steady state temperature (t1 and t2) at two points (mid point and the free end) of the one dimensional fin correspond to the minimum of the function. +// f(t1,t2) = 0.6382*t(1)^2 + 0.3191*t(2)^2 - 0.2809*t(1)*t(2) - 67.906*t(1) - 14.29*t(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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== + +clc; +// Objective function +function f = ObjectiveFunction(t) + f = 0.6382*t(1)^2 + 0.3191*t(2)^2 - 0.2809*t(1)*t(2) - 67.906*t(1) - 14.29*t(2); +endfunction +// Initial guess +x0 = [100 200]; +disp(x0, "Initial guess given to the solver is ") +input("Press enter to proceed: ") +[xopt,fopt,exitflag,output,gradient,hessian] = fminunc(ObjectiveFunction,x0) +// Result representation +clc +select exitflag +case 0 + disp("Optimal Solution Found") + disp(xopt', "The steady state temperature at the points are") + disp(fopt, "The optimum objective function value is") +case 1 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal.") + disp(xopt', "The temperature at the points are") + disp(fopt, "The objective function value is") +case 2 + disp("Maximum CPU Time exceeded. Output may not be optimal.") + disp(xopt', "The temperature at the points are") + disp(fopt, "The objective function value is") +case 3 + disp("Stop at Tiny Step.") + disp(xopt', "The temperature at the points are") + disp(fopt, "The objective function value is") +case 4 + disp("Solved To Acceptable Level.") + disp(xopt', "The temperature at the points are") + disp(fopt, "The objective function value is") +case 5 + disp("Converged to a point of local infeasibility.") +end + +disp(output) diff --git a/code/fminunc/ThreeHump-backCamel.sci b/code/fminunc/ThreeHump-backCamel.sci new file mode 100644 index 0000000..f1f2b38 --- /dev/null +++ b/code/fminunc/ThreeHump-backCamel.sci @@ -0,0 +1,69 @@ +// Reference: Urmila Diwekar, "Introduction to Aplied Optimization", 2nd Ed, Springer Science + Business Media,2008, Chapter 3 + +// Three hump-back camel function : An unconstrained problem +// F(X) = 2*X(1)^2 - 1.05*X(1)^4 + (1/6)*X(1)^6 + X(1)*X(2) + X(2)^2 +// Dimension of the problem: 2 +// The global optima for the function: F*(X) = 0 and X* = [0 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== + +clc; +// Objective function +function f = ObjectiveFunction(X) + f = 2*X(1)^2 - 1.05*X(1)^4 + (1/6)*X(1)^6 + X(1)*X(2) + X(2)^2; +endfunction +// Gradient of the function +function gf = GradientFunction(X) + gf = [4*X(1) - 4.2*X(1)^3 + X(1)^5 + X(2), X(1) + 2*X(2)]; +endfunction +// Hessian matrix of the function +function hf = HessianFunction(X) + hf = [4 - 12.6*X(1)^2 + 5*X(1)^4, 1; 1 2] +endfunction + +// Initial guess +x0 = [0,-1]; +disp(x0, "The initial guess given to the solver is") +input("Press enter to proceed: ") +// User specified parameter values +options=list("gradobj", GradientFunction,"hessian",HessianFunction); +// Calling the solver +[xopt,fopt,exitflag,output,gradient,hessian] = fminunc(ObjectiveFunction,x0,options) + +// Result representation +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/CatenaryCable_NLP.sce b/code/fsolve/CatenaryCable_NLP.sce new file mode 100644 index 0000000..1d44e94 --- /dev/null +++ b/code/fsolve/CatenaryCable_NLP.sce @@ -0,0 +1,73 @@ +// This is an example problem for solving a non linear equation + +//Reference: Steven C. Chapra. 2006. Applied Numerical Methods with MATLAB for Engineers and Scientists. McGraw-Hill Science/Engineering/Math,Chapter 6 + +//A catenary cable is one which is hung between two points not in the same vertical line. As depicted in Fig. P6.17a, it is subject to no loads other than its own weight. Thus, its weight acts as a uniformload per unit length along the cable w (N/m). A free-body diagram of a section AB is depicted in Fig. P6.17b, where TA and TB are the tension forces at the end. Based on horizontal and vertical force balances, the following differential equation model of the cable can be derived in which calculus can be employed to solve the equation for the height of the cable y as a function of distance x: +//y = (Ta/w)*cosh((w*x)/Ta) - y0 + (Ta/w) +//(a) Use a numerical method to calculate a value for the parameter TA given values for the parameters w = 10 and y0 = 5, such that the cable has a height of y = 15 at x = 50. +//(b) Develop a plot of y versus x for x = −50 to 100. + +//===================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== +clc; + +function F = CatenaryCableEqn(Ta,x,w,y0,y) +F = y - (Ta/w)*cosh((w.*x)/Ta) - y0 + (Ta/w); +endfunction + +function df = deriv_CatenaryCableEqn(Ta) + df = numderivative(CatenaryCableEqn,Ta) +endfunction + +w = 10; +y0 = 5; +y = 15; +x = 50; +Ta = 300; +Tol = 1.e-16; +[Ta_sol,v,info] = fsolve(Ta,list(CatenaryCableEqn,x,w,y0,y),deriv_CatenaryCableEqn,Tol); + +select info +case 0 + mprintf('improper input parameters'); +case 1 + mprintf('Algorithm estimates that the relative error between x and the solution is at most tol'); +case 2 + mprintf('Number of calls to fcn reached'); +case 3 + mprintf('tol is too small. No further improvement in the approximate solution x is possible'); +else + mprintf('Iteration is not making good progress'); +end + +disp(Ta_sol,'Obtained solution is ') + +input('Press enter to identify the effect of initial guess values on the optimal solution ') +Ta = -200:60:400; +nTa = length(Ta); + +for i = 1:nTa + [CurrentTa_sol(i),v,info] = fsolve(Ta(i),list(CatenaryCableEqn,x,w,y0,y),deriv_CatenaryCableEqn,Tol); + +end +disp([Ta' CurrentTa_sol],'The solution corresponding to different initial guess values') + + +input("Press enter to proceed ") +mprintf('The variation of y as the the values of x varies from -50 to 100 is plotted \n') +x = 10:10:100; +y = (Ta_sol/w)*cosh((w.*x)/Ta_sol) + y0 - (Ta_sol/w) +plot(x,y,'r*') +xlabel("x", "fontsize", 5); +ylabel("y", "fontsize", 5) + + diff --git a/code/fsolve/SystemOfNLEquations.sce b/code/fsolve/SystemOfNLEquations.sce new file mode 100644 index 0000000..f2db71f --- /dev/null +++ b/code/fsolve/SystemOfNLEquations.sce @@ -0,0 +1,57 @@ +//Reference: Edger, Hemmelblau and Lasdon,"OPtimization of chemical process", McGraw Hill, Chemical Engineering Series, 2nd edition, 2001, Appendix A + +//Solve the two non-linear equations +// X(1)^2 + X(2)^2 = 8 +// X(1)*X(2) = 4 +//===================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== + +clc; +// Set of non-linear equations +function f = NLEquations(X) + f(1) = X(1)^2 + X(2)^2 - 8 + f(2) = X(1)*X(2) - 4 +endfunction +// Multiple initial guesses. User can give their choice of values here +Xinitials = [-5 0; 5 3] +Ninitial = size(Xinitials,1); +// Precision tolerance value +Tol = 0; + +// Solving the equations with different initial guesses +for n = 1:Ninitial + Xi = Xinitials(n,:); + mprintf('Initial guess No. %d \n',n) + disp(Xi) + [X,v,info] = fsolve(Xi,NLEquations,[],Tol); + // Result representation + select info + case 0 + disp("improper input parameters.") + case 1 + disp("algorithm estimates that the relative error between x and the solution is at most tol.") + disp(X,"solution reported") + disp(v, "Value of the function at the solution") + case 2 + disp("number of calls to fcn reached") + disp(X,"solution reported") + disp(v, "Value of the function at the solution") + case 3 + disp("tol is too small. No further improvement in the approximate solution x is possible.") + disp(X,"solution obtained") + disp(v, "Value of the function at the solution") + case 4 + disp("iteration is not making good progress.") + disp(X,"solution obtained") + disp(v, "Value of the function at the solution") + end +end diff --git a/code/intfmincon/BinaryVariableProblem.sce b/code/intfmincon/BinaryVariableProblem.sce new file mode 100644 index 0000000..aa1e73a --- /dev/null +++ b/code/intfmincon/BinaryVariableProblem.sce @@ -0,0 +1,72 @@ +// M.F. Cardoso,R.L.Salcedo,S.F.D. Azevedo,D. Barbosa, A simulated anealing approach to the solution of MINLP problems, Computer and Chemical Engineering,21(12),1349-1364,1997 +// +// max f = r1*r2*r3; +// r1 = 1 - 0.1^y1*0.2^y2*0.15^y3 +// r2 = 1 - 0.05^y4*0.2^y5*0.15^y6; +// r3 = 1 - 0.02^y7*0.06^y8; +// subject to +// y1 + y2 + y3 >= 1; +// y4 + y5 + y6 >= 1; +// y7 + y8 >= 1; +// 3y1 + y2 + 2y3 + 3y4 + 2y5 + y6 + 3y7 + 2y8 <= 10; +// 0 <= y <= 1 +// All the variables (y) are integers +// Global optima: y* = [0 1 1 1 0 1 1 0]; f* = 0.94347; + +//====================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//====================================================================== + +clc; + +// Objective function +function f = ObjectiveFunction(y) + r(1) = 1 - 0.1^y(1)*0.2^y(2)*0.15^y(3); + r(2) = 1 - 0.05^y(4)*0.2^y(5)*0.15^y(6); + r(3) = 1 - 0.02^y(7)*0.06^y(8); + f = -1*prod(r); +endfunction +// Linear inequality constraints +A = [-1 -1 -1 0 0 0 0 0; 0 0 0 -1 -1 -1 0 0; 0 0 0 0 0 0 -1 -1;3 1 2 3 2 1 3 2]; +b = [-1 -1 -1 10]'; + +// Dimension of the problem +nVar = 8; +// Bounds of the problems +lb = zeros(1,8); +ub = ones(1,8); + +// Initial guess given to the solver +x0 = lb; +// indices of the integer decision variables +int = 1:nVar; +// Calling the solver +[xopt,fopt,exitflag,output,lambda] = intfmincon(ObjectiveFunction,x0,int,A,b,[],[],lb,ub) + +// Result representation +clc; +disp(x0,"Initial guess given to the solver") +select exitflag +case 0 + disp("Optimal Solution Found") + disp(xopt',"The optimum solution obtained") + disp(fopt,"The optimum value of the objective function") +case 1 + disp("Infeasible solution") +case 2 + disp("Objective Function is Continuous Unbounded") +case 3 + disp("Limit Exceeded") +case 4 + disp("User Interupt") +case 5 + disp("MINLP Errors") +end diff --git a/code/intfmincon/Pressure vessel design.sce b/code/intfmincon/Pressure vessel design.sce new file mode 100644 index 0000000..2fe9d00 --- /dev/null +++ b/code/intfmincon/Pressure vessel design.sce @@ -0,0 +1,70 @@ +// Reference: E. Sandgren,Nonlinear Integer and Discrete Programming in Mechanical I Design Optimization,Journal of Mechanical Design,JUNE 1990, Vol. 112/223 + +//A cylindrical pressure vessel is capped at both ends by hemispherical heads.The total cost, +//including the cost of material, cost of forming and welding,is to be minimized. The design variables are Ts and Th are the thicknesses of the shell and head, and R and L, the inner radius and length of the cylindrical section.These variables are denoted by X1 x2 , x3, and x4, respectively, and units for each are inches. The variables are such that R and L are continuous while Ts and Th are integer multiples of 0.0625 inch, the available thicknesses of rolled steel plates. +// Min f = 0.6224*X(1)*X(3)*X(4) + 1.7781*X(2)*X(3)^2 + 3.1661*X(1)^2*X(4) + 19.84*X(1)^2*X(3) +// subject to +// g1(X) - X1 + 0.0193 X3 < = 0 +// g2(X) - x2 + 0.00954X3 <= 0 +// g3(X) -%pi*X3^2*X4 - (4/3)*%pi*X3^3 + 1296000 <= 0 +// g4(X) X4 - 240 <= 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//====================================================================== +clc; + +// Onjective fucntion +function f = ObjectiveFunction (X) + X(1:2) = X(1:2)*0.0625; + f = 0.6224*X(1)*X(3)*X(4) + 1.7781*X(2)*X(3)^2 + 3.1661*X(1)^2*X(4) + 19.84*X(1)^2*X(3); +endfunction +// Non linear equality and inequality constraints +function [C,Ceq] = NLconstraints(X) + X(1:2) = (X(1:2))*0.0625; + C = -%pi*X(3)^2*X(4) - (4/3)*%pi*X(3)^3 + 1296000; + Ceq = []; +endfunction +// Linear inequality constraints +A = [-0.0625 0 0.0193 0;0 -0.0625 0.00954 0;0 0 0 1]; +b = [0 0 240]'; +// Bounds of the variables +lb = [1 1 10 10]; +ub = [99 99 200 200]; +nVar = length(lb); +// Initial guess given to the solver +x0 = [20 10 58.291 43.69]; +// indices of the integer decision variables +int = [1 2]; +// Calling the solver +[xopt,fopt,exitflag,output,lambda] = intfmincon(ObjectiveFunction,x0,int,A,b,[],[],lb,ub,NLconstraints) + +// Result representation +// Converting the integer variables to the discrete variable +x0(1:2) = x0(1:2)*0.0625; +clc; +disp(x0,"Initial guess given to the solver") +select exitflag +case 0 + disp("Optimal Solution Found") + disp(xopt',"The optimum solution obtained") + disp(fopt,"The optimum value of the objective function") +case 1 + disp("Converged to a point of local infeasibility") +case 2 + disp("Objective Function is Continuous Unbounded") +case 3 + disp("Limit Exceeded") +case 4 + disp("User Interupt") +case 5 + disp("MINLP Errors") +end diff --git a/code/intqpipopt/intqpipoptExample1.sce b/code/intqpipopt/intqpipoptExample1.sce new file mode 100644 index 0000000..cb5c648 --- /dev/null +++ b/code/intqpipopt/intqpipoptExample1.sce @@ -0,0 +1,76 @@ +//Reference: K. Deep et al.,"A real coded genetic algorithm for solving integer and mixed integer optimization problems", Applied Mathematics and Computation, 212, p 505-518,2009 +// min f(x) = x1^2 + x2^2 +x3^2 +x4^2 +x5^2 -8x1 - 2x2 - 3x3 - x4 -2x5 +// Subject to +// x1 + x2 + x3 + x4 + x5 <= 400; +// x1 + 2x2 + 2x3 + x4 + 6x5 <= 800 +// 2x1 + x2 + 6x3 <= 200 +// x3 + x4 + 5x5 <= 200 +// x1 + x2 + x3 + x4 + x5 >= 55 +// x1 + x2 + x3 + x4 >= 48 +// x2 + x4 + x5 >= 34 +// 6x1 + 7x5 >= 104 +// The known optimal solution is f*(x) = 807 at x* = [16. 22. 5. 5. 7.] +//===================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== + +clc; +// Number of decision variables +nVar = 5; +// Num\ber of constraints +nCon = 8; +//integer constraints +intcon = 1:5; +//coefficients of the linear term +f = [-8 -2 -3 -1 -2]; +// Hessian matrix +H = [2 0 0 0 0; +0 2 0 0 0 +0 0 6 0 0 +0 0 0 8 0 +0 0 0 0 4]; +// Bounds of the problem +lb = zeros(1,nVar); +ub = 99*ones(1,nVar); +// Constraint matrix +A = [1 1 1 1 1; +1 2 2 1 6; +2 1 6 0 0; +0 0 1 1 5; +-1 -1 -1 -1 -1 +-1 -1 -1 -1 0; +0 -1 0 -1 -1; +-6 0 0 0 -7]; +b = [400 800 200 200 -55 -48 -34 -104]'; +// Initial guess to the solver +x0 = lb; +// Calling the solver +[xopt,fopt,exitflag,output] = intqpipopt(H,f,intcon,A,b,[],[],lb,ub,x0) + +// Result representation +disp(x0, "The initial guess given to the solver") + +select exitflag +case 0 + disp("Optimal Solution Found") + disp(xopt',"The optimal solution determined by the solver") + disp(fopt,"The optimal objective function") +case 1 + disp("InFeasible Solution") +case 2 + disp("Objective Function is Continuous Unbounded") +case 3 + disp("Limit Exceeded") +case 4 + disp("User Interrupt") +case 5 + disp("MINLP Error") +end diff --git a/code/intqpipopt/intqpipoptExample2.sce b/code/intqpipopt/intqpipoptExample2.sce new file mode 100644 index 0000000..c10e049 --- /dev/null +++ b/code/intqpipopt/intqpipoptExample2.sce @@ -0,0 +1,63 @@ +//Reference: K. Deep et al.,"A real coded genetic algorithm for solving integer and mixed integer optimization problems", Applied Mathematics and Computation, 212, p 505-518,2009 +//===================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//===================================================================== + +clc; +// Number of decision variables +nVar = 5; +// Num\ber of constraints +nCon = 6; +//integer constraints +intcon = 1:5; +//coefficients of the linear term +f = zeros(1,nVar); +// Hessian matrix +H = [2 0 0 0 0; +0 2 0 0 0 +0 0 2 0 0 +0 0 0 2 0 +0 0 0 0 2]; +// Bounds of the problem +lb = zeros(1,nVar); +ub = 3*ones(1,nVar); +// Constraint matrix +A = [-1 -2 0 -1 0; +0 -1 -2 0 0; +-1 0 0 0 -2; +1 2 2 0 0; +2 0 1 0 0; +1 0 0 0 4]; +b = [-4 -3 -5 6 4 13]'; +// Initial guess to the solver +x0 = lb; +// Calling the solver +[xopt,fopt,exitflag,output] = intqpipopt(H,f,intcon,A,b,[],[],lb,ub,x0) + +// Result representation +disp(x0, "The initial guess given to the solver") + +select exitflag +case 0 + disp("Optimal Solution Found") + disp(xopt',"The optimal solution determined by the solver") + disp(fopt,"The optimal objective function") +case 1 + disp("InFeasible Solution") +case 2 + disp("Objective Function is Continuous Unbounded") +case 3 + disp("Limit Exceeded") +case 4 + disp("User Interrupt") +case 5 + disp("MINLP Error") +end diff --git a/code/linprog/AirforceBombDuality.sce b/code/linprog/AirforceBombDuality.sce new file mode 100644 index 0000000..3bc2ef3 --- /dev/null +++ b/code/linprog/AirforceBombDuality.sce @@ -0,0 +1,73 @@ +//Reference: J.K. Sharma,"OPerations Research Theory and Applications", Macmillan Publishers India Ltd., 5th edition, 2013, Chapter 4. + +//An Air Force is experimenting with three types of bombs P,Q and R in which three kinds of explosives, viz., A, B and C will be used. Taking the various factors into account, it has been decided to use the maximum 600 kg of explosive A, at least 480 kg of B and exactly 540 kg of explosive C. Bomb P requires 3,2,2 kg, bomb Q requires 1,4,3 kg and bomb R requires 4,2,3 kg of explosives A,B and C respectively. Bomb P is estimated to give the equivalent of a 2 ton explosion, bomb Q, a 3 ton of explosion and bomb R, a 4 ton of explosion respectively. Under what production schedule can the Air Force make the biggest bang? +//============================================================================= +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//============================================================================ + +clc; +// =============================Primal Model=================================== +nVar = 3; + +// Linear inequalities +A = [3 1 4; -2 -4 -2]; +b = [600 -480]'; +nIneqCon = length(b); + +// Linear equality +Aeq = [2 3 3]; +beq = 540; +nEqCon = length(beq); + +lb = zeros(1,nVar); +ub = []; + +Cost = -[2 3 4 ]'; + +// =============================Dual Model=================================== + +// Inequality constraints +dual_A = -[A;Aeq]'; +dual_beq = Cost; +n_dualIneq = length(dual_beq); +Coeff_Slack = [eye(n_dualIneq,n_dualIneq);]; +dual_Aeq = [dual_A Coeff_Slack]; +dual_Cost = [b;beq; zeros(n_dualIneq,1)]; + +dual_nVar = length(dual_Cost); +lb = [zeros(1,nIneqCon) -%inf*ones(nEqCon) zeros(1,(dual_nVar-nIneqCon-nEqCon))]; +ub = []; + +[xopt,fopt,exitflag,output,lambda]=linprog(dual_Cost, [],[], dual_Aeq, dual_beq, lb,ub); +clc + +select exitflag +case 0 + disp(xopt(1:nIneqCon+nEqCon)',"Dual values of the primal constraints") + disp(xopt(nIneqCon+nEqCon+1:dual_nVar)',"Dual values of the primal variables") + TopRow = ["Variable " "Rate of change in optima"]; + Table = string([(1:length(xopt))' -xopt]); + disp("The change in optima due to one unit increase in each variables are given as below") + disp([TopRow;Table]) + disp(["The optimal cost is ", string(fopt)]) +case 1 + disp("Primal Infeasible") +case 2 + disp("Dual Infeasible") +case 3 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 4 + disp("Solution Abandoned") +case 5 + disp("Primal objective limit reached") +case 6 + disp("Dual objective limit reached") +end diff --git a/code/linprog/MoulderProblemDuality.sce b/code/linprog/MoulderProblemDuality.sce new file mode 100644 index 0000000..c7a13c3 --- /dev/null +++ b/code/linprog/MoulderProblemDuality.sce @@ -0,0 +1,87 @@ +// Example problem for sensitivity analysis - using duality theorem to determine the reduced cost + +//Reference : Bradley, Hax, and Magnanti,"Applied Mathematical Programming", Addison-Wesley, 1977, chapter 3 + +//A custom molder has one injection-molding machine and two different dies to fit the machine. Due to differences in number of cavities and cycle times, withthe first die he can produce 100 cases of six-ounce juice glasses in six hours,while with the second die he can produce 100 cases of ten-ounce fancy cocktail glasses in five hours. The molder is approached by a new customer to produce a champagne glass. The production time for the champagne glass is 8 hours per hundred cases. He prefers to operate only on a schedule of 60 hours of production per week. He stores the week’s production in his own stockroom where he has an effective capacity of 15,000 cubic feet. A case of six-ounce juice glasses requires 10 cubic feet of storage space, while a case of ten-ounce cocktail glasses requires 20 cubic feet due to special packaging.The storage space required for the champagne glasses is 10 cubic feet per one case. The contribution of the six-ounce juice glasses is $5.00 per case; however, the only customer available will not accept more than 800 cases per week. The contribution of the ten-ounce cocktail glasses is $4.50 per case and there is no limit on the amount that can be sold where as the contribution of champagne glasses is $6.00 per case with no limit on the bound. How many cases of each type of glass should be produced each week in order to maximize the total contribution? +//============================================================================= +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//============================================================================ + +clc; +// -----------------------------Primal model------------------------------- +nProducts = 3; +TotalProcessingTime = 60; +RequiredTime = [6 5 8]; +StokeRoomCapacity = 15000; +RequiredSpace = [10 20 10]*100; +Demand = 8; +Revenue = [5 4.5 6]; +nVar = 3 // Each variable indicates number of hundreds of units of each product + +// Linear constraints +// Processing time constraint +A1 = RequiredTime; +b1 = TotalProcessingTime; +// Space constraint +A2 = RequiredSpace; +b2 = StokeRoomCapacity; +// Demand constraint +nDemand = length(Demand); +A3 = eye(nDemand,nVar); +b3 = Demand'; + +A = [A1;A2;A3]; b = [b1;b2;b3]; + +//Addition of slack variables +nInEqConstraints = length(b); +Aeq = eye(nInEqConstraints,nInEqConstraints); +Aeq = [A Aeq]; +beq = b; + +Cost = -[Revenue zeros(1,nInEqConstraints)]'; + +// -------------------------------Dual model------------------------------ +dualA = -A'; +dualb = Cost(1:nVar); + +dual_nInEqConstraints = length(dualb); +dualAeq = eye(dual_nInEqConstraints,dual_nInEqConstraints); +dualAeq = [dualA dualAeq]; +dualbeq = dualb; +dualCost = beq; +dualCost = [dualCost; zeros(dual_nInEqConstraints,1)]; +dual_nVar = length(dualCost); +dual_lb = zeros(1,size(dualAeq,2)); + +[dual_xopt,dual_fopt,dual_exitflag,dual_output,dual_lambda] = linprog(dualCost,[],[], dualAeq, dualbeq, dual_lb,[]); +clc; +select dual_exitflag +case 0 + disp(dual_xopt(1:nInEqConstraints)',"Dual values of the primal constraints") + disp(dual_xopt(nInEqConstraints+1:dual_nVar)',"Dual values of the primal variables") + TopRow = ["Variable " "Rate of change in optima"]; + Table = string([(1:length(dual_xopt))' -dual_xopt]); + disp("The change in optima due to one unit increase in each variables are given as below") + disp([TopRow;Table]) + disp(["The optimal cost is ", string(dual_fopt)]) +case 1 + disp("Primal Infeasible") +case 2 + disp("Dual Infeasible") +case 3 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 4 + disp("Solution Abandoned") +case 5 + disp("Primal objective limit reached") +case 6 + disp("Dual objective limit reached") +end diff --git a/code/linprog/MoulderProblem_RHSchange.sce b/code/linprog/MoulderProblem_RHSchange.sce new file mode 100644 index 0000000..95721c5 --- /dev/null +++ b/code/linprog/MoulderProblem_RHSchange.sce @@ -0,0 +1,78 @@ +// Example problem of sensitivity analysis - changes in RHS + +//Reference : Bradley, Hax, and Magnanti,"Applied Mathematical Programming", Addison-Wesley, 1977, chapter 3 + +//A custom molder has one injection-molding machine and two different dies to fit the machine. Due to differences in number of cavities and cycle times, withthe first die he can produce 100 cases of six-ounce juice glasses in six hours,while with the second die he can produce 100 cases of ten-ounce fancy cocktail glasses in five hours. The molder is approached by a new customer to produce a champagne glass. The production time for the champagne glass is 8 hours per hundred cases. He prefers to operate only on a schedule of 60 hours of production per week. He stores the week’s production in his own stockroom where he has an effective capacity of 15,000 cubic feet. A case of six-ounce juice glasses requires 10 cubic feet of storage space, while a case of ten-ounce cocktail glasses requires 20 cubic feet due to special packaging.The storage space required for the champagne glasses is 10 cubic feet per one case. The contribution of the six-ounce juice glasses is $5.00 per case; however, the only customer available will not accept more than 800 cases per week. The contribution of the ten-ounce cocktail glasses is $4.50 per case and there is no limit on the amount that can be sold where as the contribution of champagne glasses is $6.00 per case with no limit on the bound. How many cases of each type of glass should be produced each week in order to maximize the total contribution? +//============================================================================= +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//============================================================================ + +clc; + +nProducts = 3; +TotalProcessingTime = 60; +RequiredTime = [6 5 8]; +StokeRoomCapacity = 15000; +RequiredSpace = [10 20 10]*100; +Demand = 8; +Revenue = [5 4.5 6]; +nVar = nProducts // Each variable indicates number of hundreds of units of each product + +// Linear constraints +// Processing time constraint +A1 = RequiredTime; +b1 = TotalProcessingTime; +// Space constraint +A2 = RequiredSpace; +b2 = StokeRoomCapacity; +// Demand constraint +nDemand = length(Demand); +A3 = eye(nDemand,nVar); +b3 = Demand'; + +A = [A1;A2;A3]; b = [b1;b2;b3]; + +//Addition of slack variables +nInEqConstraints = length(b); +Aeq = eye(nInEqConstraints,nInEqConstraints); +Aeq = [A Aeq]; +beq = b; + +Cost = -[Revenue zeros(1,nInEqConstraints)]'; +lb = zeros(1,size(Aeq,2)); + +[xopt,fopt,exitflag,output,lambda]=linprog(Cost, [],[], Aeq, beq, lb,[]); +clc; +select exitflag +case 0 + mprintf('The optimum solution has been found \n') + disp(xopt(1:nVar)', "The optimum number of six-ounce juice glasses,ten-ounce fancy cocktail glasses and champagne glass") + SlackVariables = xopt(nVar+1:length(xopt)); + disp(SlackVariables',"The surplus resource available") + Index = find(SlackVariables>0); + newbeq = beq(Index) - SlackVariables(Index); + changedbeq = beq; + changedbeq(Index) = newbeq; + mprintf(' The value beyond %f of the RHS value of constraint number %d will not change the optima. \n ',newbeq,Index) + mprintf('\n The optimal cost is %f',fopt) +case 1 + disp("Primal Infeasible") +case 2 + disp("Dual Infeasible") +case 3 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 4 + disp("Solution Abandoned") +case 5 + disp("Primal objective limit reached") +case 6 + disp("Dual objective limit reached") +end diff --git a/code/linprog/MultiStage_Planning_Problem.sce b/code/linprog/MultiStage_Planning_Problem.sce new file mode 100644 index 0000000..51d37dd --- /dev/null +++ b/code/linprog/MultiStage_Planning_Problem.sce @@ -0,0 +1,165 @@ +//Reference : Bradley, Hax, and Magnanti,"Applied Mathematical Programming", Addison-Wesley, 1977, chapter 5 + +// The data given in the Multi-stage planning problem are as follows +//An automobile tire company has the ability to produce both nylon and fiberglass tires. During the next three months they have agreed to deliver tires as follows +// ========================================= +// Date Nylon Fiberglass +// ----------------------------------------- +// June 30 4000 1000 +// July 31 8000 5000 +// August 31 3000 5000 +// ========================================= +//The company has two presses, a Wheeling machine and a Regal machine, and appropriate molds that can be used to produce these tires, with the following production hours available in the upcoming months +// ================================================ +// Wheeling Machine Regal Machine +// ------------------------------------------------ +// June 700 1500 +// July 300 400 +// August 1000 300 +// ================================================ +//The production rates for each machine-and-tire combination, in terms of hours per tire, are as follows: +//====================================================== +// Wheeling Machine Regal Machine +//------------------------------------------------------ +//Nylon 0.15 0.16 +//Fiberglass 0.12 0.140 +//====================================================== +//The variable costs of producing tires are $5.00 per operating hour, regardlessof which machine is being used or which tire is being produced. There is also an inventory-carrying charge of $0.10 per tire per month. How should the production be scheduled in order to meet the delivery requirements at minimum costs? +//Reported solution: +//X = [1867 7633 3500 0 5500 2500 0 2500 2500 0 0 2667 333 5000 0]; +//=========================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//============================================================================ +clc; +nProducts = 2; +nMachines = 2; +nPeriods = 3; +Demand = [4000 1000;8000 5000; 3000 5000]; +AvailMachineTime = [700 1500;300 400; 1000 300]; +ProdRate = [0.15 0.16;0.12 0.14]; +OperatingCost = 5; +InventoryCost = 0.1; + +nVar = nProducts*nMachines*nPeriods + (nPeriods-1)*nProducts; // Dimension of the problem is determined + +// Linear equality constraints +// Demand constraints +nEqConstraints = nProducts*nPeriods +Aeq = zeros(nEqConstraints,nVar); +// Demand constraints for period 1 +Aeq1 = zeros(nProducts,nVar); +for i = 1:nProducts + index1 = (i-1)*nProducts+1:i*nProducts; + Aeq1(i,index1) = 1; + index2 = nMachines*nProducts+i; + Aeq1(i,index2) = -1; + beq1(i,1) = Demand(1,i); +end +// Demand constraints for period 2 to (nPeriods-1)th period +Aeq2 = zeros(nProducts,nVar); +for i = 2:nPeriods-1 + for j = 1:nProducts + index3 = (i-1)*(nProducts*nMachines+nProducts)+(j-1)*nMachines+1:(i-1)*(nProducts*nMachines+nProducts)+(j-1)*nMachines+nMachines; + Aeq2(j,index3) = 1; + index4 = (i-1)*(nProducts*nMachines+nProducts) - nProducts+j; + Aeq2(j,index4) = 1; + index5 = i*(nProducts*nMachines+nProducts)- nProducts+j + Aeq2(j,index5) = -1; + beq2(j,1) = Demand(i,j); + end +end + +// Demand constraints for last period +Aeq3 = zeros(nProducts,nVar); +for i = 1:nProducts + index6 = (nProducts*nMachines+nProducts)*(nPeriods-1)+(i-1)*nProducts+1:(nProducts*nMachines+nProducts)*(nPeriods-1)+i*nProducts + Aeq3(i,index6) = 1; + index7 = (nProducts*nMachines+nProducts)*(nPeriods-1) - nProducts+i; + Aeq3(i,index7) = 1; + beq3(i,1) = Demand(nPeriods,i); +end +Aeq = [Aeq1;Aeq2;Aeq3]; +beq = [beq1;beq2;beq3]; + +// Linear inequality constraints +// Machine time constraints +for i = 1:nPeriods + for j = 1:nProducts + Cindex = (i-1)*(nProducts*nMachines+nProducts)+j:nProducts:(i-1)*(nProducts*nMachines+nProducts)+nMachines+j; + Rindex = (i-1)*nProducts+j; + A(Rindex,Cindex) = ProdRate(:,j)'; + b(Rindex,1) = AvailMachineTime(i,j); + end +end + +// Objective function +TotalProductionCost = []; +for j = 1:nProducts + TotalProductionCost = [TotalProductionCost ProdRate(j,:)*OperatingCost]; +end + +for i = 1:nPeriods + index = (i-1)*(nProducts*nMachines+nProducts)+1:(i-1)*(nProducts*nMachines+nProducts)+nProducts*nMachines; + nindex = length(index); + cost(index,1) = TotalProductionCost'; + cost(index(nindex)+1:index(nindex)+nProducts,1) = InventoryCost; +end +cost(nVar+1:nVar+nProducts,1) = []; +lb = zeros(1,nVar); + +[xopt,fopt,exitflag,output,lambda]=linprog(cost, A, b, Aeq, beq, lb,[]); + +//Result representation + +select exitflag +case 0 + disp(" Optimal Solution Found") + M = [" "]; + for m = 1:nMachines + M = [M strcat(["Machine",string(m)])]; + end + P = []; + for p = 1:nProducts + P = [P;strcat(["Product ",string(p)])]; + end + + for i = 1:nPeriods + Sol = []; + for j = 1:nProducts + Ind1 = (i-1)*(nProducts*nMachines + nProducts)+(j-1)*nProducts+1:(i-1)*(nProducts*nMachines + nProducts)+j*nProducts; + Sol = [Sol;xopt(Ind1)']; + end + + disp(strcat(["Production schedule for the Period ", string(i)])); + disp([M; [P string(Sol)]]); + end + + for i = 1:nPeriods-1 + ind = i*(nProducts*nMachines+1:nProducts*nMachines+nProducts); + inventory = xopt(ind); + disp(strcat(["Inventory at the Period ", string(i)])); + disp([P string(inventory)]); + end + + disp(["The optimal cost is ", string(fopt)]) +case 1 + disp("Primal Infeasible") +case 2 + disp("Dual Infeasible") +case 3 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 4 + disp("Solution Abandoned") +case 5 + disp("Primal objective limit reached") +case 6 + disp("Dual objective limit reached") +end diff --git a/code/linprog/ProductMix.sce b/code/linprog/ProductMix.sce new file mode 100644 index 0000000..5d66fed --- /dev/null +++ b/code/linprog/ProductMix.sce @@ -0,0 +1,87 @@ +//Refernce: H.P. Williams,"Model building in mathematical programming",Fifth Edition,AJohn Wiley & Sons, Ltd., Publication,2013, chapter 1.2 + +//An engineering factory can produce five types of product (PROD 1, PROD 2,... , PROD 5) by using two production processes: grinding and drilling. After deducting raw material costs, each unit of each product yields the following contributions to profit: +// ===================================================== +// PROD1 PROD2 PROD3 PROD4 PROD5 +// 550 600 350 400 200 +// ===================================================== +//Each unit requires a certain time on each process. These are given below (in +//hours). A dash indicates when a process is not needed. +// ================================================================= +// PROD1 PROD2 PROD3 PROD4 PROD5 +// ----------------------------------------------------------------- +// Grinding 12 20 - 25 15 +// Drilling 10 8 16 - - +// ================================================================= +//In addition, the final assembly of each unit of each product uses 20 hours of +//an employee’s time. The factory has three grinding machines and two drilling machines and works a six-day week with two shifts of 8 hours on each day. Eight workers are employed in assembly, each working one shift a day. The problem is tofind how much of each product is to be manufactured so as to maximize the totalprofit contribution. The minimum demand for each product is added to this problem and details are as follows +// ================================================================= +// PROD1 PROD2 PROD3 PROD4 PROD5 +// ----------------------------------------------------------------- +// Demand 5 5 2 7 2 +// ================================================================= + +//============================================================================= +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//============================================================================ + +clc; +nProducts = 5; +nMachines = [3 2]; +nWorkers = 8; +nDays = 6; +nHours = 8; +nMachineShift = 2; +nWorkerShift = 1; +prodProfit = [550 600 350 400 200]'; +ProcessingTime = [12 20 0 25 15; 10 8 16 0 0; 20*ones(1,nProducts)]; + +// Linear inequalities +A1 = ProcessingTime; +b1 = (nDays*nHours)*[nMachines*nMachineShift nWorkers*nWorkerShift]'; +// Demand constraints +A2 = [-1 0 0 0 0;0 -1 0 0 0;0 0 -1 0 0; 0 0 0 -1 0; 0 0 0 0 -1]; +b2 = [-5 -5 -2 -7 -2]'; + +A = [A1;A2]; +b = [b1;b2]; +lb = zeros(1,nProducts); +Cost = -prodProfit; +[xopt,fopt,exitflag,output,lambda]=linprog(Cost, A, b, [], [], lb,[]); + +// Result representation +select exitflag +case 0 + disp(" Optimal Solution Found") + P = []; + for p = 1:nProducts + P = [P strcat(["PROD",string(p)])]; + end + + disp([P;string(xopt')]); + disp(["The optimal cost is ", string(fopt)]) +case 1 + disp("Primal Infeasible") + nViolatedConstraints = sum(A*xopt-b>0) + if nViolatedConstraints + disp(strcat(["No of constraints violated : ", string(nViolatedConstraints)])); + end +case 2 + disp("Dual Infeasible") +case 3 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 4 + disp("Solution Abandoned") +case 5 + disp("Primal objective limit reached") +case 6 + disp("Dual objective limit reached") +end diff --git a/code/linprog/Transportation_Problem.sce b/code/linprog/Transportation_Problem.sce new file mode 100644 index 0000000..9df7c49 --- /dev/null +++ b/code/linprog/Transportation_Problem.sce @@ -0,0 +1,121 @@ +//Reference : Dantzig G. B., Linear Programming and Extensions, Princeton University Press, Princeton, New Jersey, 1963, Chapter 3-3. + +//The objective of the problem is to minimize the transportation cost of goods from the canneries to the warehouses,by satisfying the supply and demand constraints. Inorder to simplify the problem, it is assumed that there are two canneries and three warehouses. The supply data of canneries are as follows +// ============================= +// Supply +// Cannery I 350 +// Cannery II 650 +// ============================= +//The demands at the warehouses are as given +// ============================= +// Demand +// Warehouse A 300 +// Warehouse B 300 +// Warehouse C 300 +// ============================= +//ShippingCost between the canneries and warehouses are +// ========================================================= +// Warehouse A Warehouse B Warehouse C +// ========================================================= +// Cannery I 2.5 1.7 1.8 +// Cannery II 2.5 1.8 1.4 +// ========================================================= +// It should be noted that, this case is solved without considering the inventory. +// +//=========================================================================== +// 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: Remya Kommadath +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +//============================================================================ + +clc; + +nCanneries = 2; +nWarehouses = 3; +ShippingCost = [2.5 1.7 1.8; 2.5 1.8 1.4]; +demand = [300 300 300]'; +supply = [350 650]'; + +nVar = nCanneries*nWarehouses; +lb = zeros(1,nVar); // lower bounds of the amount of commodity to be shipped from each canneries to warehouses; + +// Linear inequality constraints +// Total supply to the warehouses from a cannery should not exceed the corresponding supply limit of the same cannery +for i = 1:nCanneries + indeX = (i-1)*nWarehouses+1:i*nWarehouses; + As(i,indeX) = ones(1,nWarehouses); +end +Bs = supply; +As +// Total supply to a warehouses from the canneries should satisfy the corresponding demand of the same warehouses +//Ad = zeros() +for i = 1:nWarehouses + indeX = i:nWarehouses:nVar; + Ad(i,indeX) = -1; +end +Bd = -demand; + +A = cat(1,As,Ad); +b = cat(1,Bs,Bd); + +// Minimize the total transportation cost between all plants and markets +Cost = zeros(nVar,1); +for i = 1:nCanneries + indeX = (i-1)*nWarehouses+1:i*nWarehouses; + Cost(indeX) = (ShippingCost(i,:)'); +end + +[xopt,fopt,exitflag,output,lambda] = linprog(Cost,A,b,[],[],lb,[]); + +// Display of results +select exitflag +case 0 + disp(" Optimal Solution Found") + // xopt is transformed into a nCanneries X nWarehouses sized string matrix + for i = 1:nCanneries + indeX = (i-1)*nWarehouses+1:i*nWarehouses; + X(i,1:nWarehouses) = xopt(indeX,1)'; + end + Values = string(X); + + // Labelling each markets as M_1,M_2,...,M_nWarehouses + Markets(1,1) = [" "]; + for j = 2:nWarehouses+1 + Markets(1,j) = strcat(["M_",string(j-1)]); + end + + // Labelling each plants as P_1,P_2,...,P_nCanneries + for i = 1:nCanneries + Plants(i) = strcat(["P_",string(i)]); + end + + table = [Markets; [Plants Values]]; + f = gcf(); + clf + as = f.axes_size; // [width height] + ut = uicontrol("style","table",.. + "string",table,.. + "position",[150 as(2)-220 300 87]) // => @top left corner of figure + + + // Modify by hand some values in the table. Then get them back from the ui: + matrix(ut.string,size(table)) +case 1 + disp("Primal Infeasible") +case 2 + disp("Dual Infeasible") +case 3 + disp("Maximum Number of Iterations Exceeded. Output may not be optimal") +case 4 + disp("Solution Abandoned") +case 5 + disp("Primal objective limit reached") +case 6 + disp("Dual objective limit reached") +end |