diff --git a/picat/diet2.pi b/picat/diet2.pi new file mode 100644 index 00000000..e3e457b7 --- /dev/null +++ b/picat/diet2.pi @@ -0,0 +1,102 @@ +/* + + Diet problem in Picat. + + From the GLPK model diet.mod + """ + STIGLER'S NUTRITION MODEL + + This model determines a least cost diet which meets the daily + allowances of nutrients for a moderately active man weighing 154 lbs. + + References: + Dantzig G B, "Linear Programming and Extensions." + Princeton University Press, Princeton, New Jersey, 1963, + Chapter 27-1. + """ + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +go => + a(A), + b(B), + NumNutr = A[1].length, + NumFood = A.length, + + X = new_list(NumFood), + X :: 0.0..1.0, + + Cost #= sum(X), + Cost #>= 0, + + % nutrient balance (units) + foreach(N in 1..NumNutr) + sum([A[F,N] * X[F] : F in 1..NumFood]) #= B[N] + end, + + + solve($[min(Cost)], X), + + printf("\nTotal cost: $%0.4f\n", Cost), + % println(x=X), + println("Dollars of food f to be purchased daily:"), + Food = ["Wheat", "Cornmeal", "Cannedmilk", "Margarine", "Cheese", "Peanut-B", "Lard", + "Liver", "Porkroast", "Salmon", "Greenbeans", "Cabbage", "Onions", "Potatoes", "Spinach", + "Sweet-Pot", "Peaches", "Prunes", "Limabeans", "Navybeans"], + + foreach({F,C} in zip(Food, X)) + printf("%-10w: $%0.4f\n", F, C) + end, + + nl. + + + +% Calorie 3 % thousands +% Protein 70 % grams +% Calcium 0.8 % grams +% Iron 12 % milligrams +% Vitamin-A 5 % thousands IUs +% Vitamin-B1 1.8 % milligrams +% Vitamin-B2 2.7 % milligrams +% Niacin 18 % milligrams +% Vitamin-C 75 % milligrams ; + +b(B) => B = [3.0, 70.0, 0.8, 12.0, 5.0, 1.8, 2.7, 18.0, 75.0]. + + +a(A) => +A = [ +% Calorie Protein Calcium Iron Vitamin-A Vitamin-B1 Vitamin-B2 Niacin Vitamin-C +% (1000) (g) (g) (mg) (1000IU) (mg) mg) (mg) (mg) +[44.7, 1411.0, 2.0, 365.0, 0.0, 55.4, 33.3, 441.0, 0.0], +[36.0, 897.0, 1.7, 99.0, 30.9, 17.4, 7.9, 106.0, 0.0], +[ 8.4, 422.0, 15.1, 9.0, 26.0, 3.0, 23.5, 11.0, 60.0], +[20.6, 17.0, 0.6, 6.0, 55.8, 0.2, 0.0, 0.0, 0.0], +[ 7.4, 448.0, 16.4, 19.0, 28.1, 0.8, 10.3, 4.0, 0.0], +[15.7, 661.0, 1.0, 48.0, 0.0, 9.6, 8.1, 471.0, 0.0], +[41.7, 0.0, 0.0, 0.0, 0.2, 0.0, 0.5, 5.0, 0.0], +[ 2.2, 333.0, 0.2, 139.0, 169.2, 6.4, 50.8, 316.0, 525.0], +[ 4.4, 249.0, 0.3, 37.0, 0.0, 18.2, 3.6, 79.0, 0.0], +[ 5.8, 705.0, 6.8, 45.0, 3.5, 1.0, 4.9, 209.0, 0.0], +[ 2.4, 138.0, 3.7, 80.0, 69.0, 4.3, 5.8, 37.0, 862.0], +[ 2.6, 125.0, 4.0, 36.0, 7.2, 9.0, 4.5, 26.0, 5369.0], +[ 5.8, 166.0, 3.8, 59.0, 16.6, 4.7, 5.9, 21.0, 1184.0], +[14.3, 336.0, 1.8, 118.0, 6.7, 29.4, 7.1, 198.0, 2522.0], +[ 1.1, 106.0, 0.0, 138.0, 918.4, 5.7, 13.8, 33.0, 2755.0], +[ 9.6, 138.0, 2.7, 54.0, 290.7, 8.4, 5.4, 83.0, 1912.0], +[ 8.5, 87.0, 1.7, 173.0, 86.8, 1.2, 4.3, 55.0, 57.0], +[12.8, 99.0, 2.5, 154.0, 85.7, 3.9, 4.3, 65.0, 257.0], +[17.4, 1055.0, 3.7, 459.0, 5.1, 26.9, 38.2, 93.0, 0.0], +[26.9, 1691.0, 11.4, 792.0, 0.0, 38.4, 24.6, 217.0, 0.0]]. + diff --git a/picat/dirichlet.pi b/picat/dirichlet.pi new file mode 100644 index 00000000..fb504491 --- /dev/null +++ b/picat/dirichlet.pi @@ -0,0 +1,52 @@ +/* + + Triangle problem in Picat. + + From the lecture notes Modelling with Constraints + http://www.cse.unsw.edu.au/~cs4418/2008/Lectures/Modelling2.ppt, page 45ff + + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + +import mip. + +main => go. + +go => + temp(Temp), + N = Temp.length, + + % Note .vars() here. + Temp.vars() :: 0.0..100.0, + + foreach(I in 2..N-1,J in 2..N-1) + % Temp[I,J] #= 1.0/4.0*(Temp[I-1,J] + Temp[I+1,J] + Temp[I,J-1] + Temp[I,J+1]) + Temp[I,J]*4.0 #= (Temp[I-1,J] + Temp[I+1,J] + Temp[I,J-1] + Temp[I,J+1]) + end, + + solve(Temp), + + foreach(Row in Temp) + foreach(T in Row) + printf("%10.3f ", T) + end, + nl + end, + + HalfN = N div 2, + printf("Tempature at Temp[%d,%d]: %0.3f\n", HalfN, HalfN, Temp[HalfN,HalfN]), + + nl. + +temp(Temp) => + Temp = + [[100.0,100.0,100.0,100.0,100.0,100.0,100.0], + [ 0.0, _, _, _, _, _, 0.0], + [ 0.0, _, _, _, _, _, 0.0], + [ 0.0, _, _, _, _, _, 0.0], + [ 0.0, _, _, _, _, _, 0.0], + [ 0.0, _, _, _, _, _, 0.0], + [ 0.0, _, _, _, _, _, 0.0]]. \ No newline at end of file diff --git a/picat/index.html b/picat/index.html index 0fe2b854..04c8fa68 100644 --- a/picat/index.html +++ b/picat/index.html @@ -36,7 +36,7 @@

My Picat programs/models

They are separated in different sections: - -

Project Euler problem

Here are some solutions for the Project Euler problems:
@@ -391,7 +401,7 @@

Project Euler problem

  • euler39.pi: Euler #39
  • euler40.pi: Euler #40
  • euler41.pi: Euler #41 -
  • euler41.pi: Euler #42 (data file words.txt) +
  • euler41.pi: Euler #42
  • euler43.pi: Euler #43
  • euler44.pi: Euler #44
  • euler45.pi: Euler #45 @@ -565,13 +575,12 @@

    Miscellaneous programs

  • tic_tac_bench.pi -
  • allpartitions.pi: Generate all partitions of a number
  • arbitrage_loops.pi: Arbitrage loops. Cf bitcoin.pi
  • analogy.pi: Evan's Analogy program (ported from Sterling "The Art of Prolog")
  • apl_util.pi: APL utils. Some utilities inspired by APL and/or K. (Is now a module in Picat distribution.)
  • base_conversion.pi: Base conversion.
  • beds1.pi: Bunk beds puzzle -
  • bitcoind.pi: Reading bitcoin currency and make an optimal arbitrage loop. Cf arbitrage_loops.pi. +
  • bitcoin.pi: Reading bitcoin currency and make an optimal arbitrage loop. Cf arbitrage_loops.pi.
  • bridge_party.pi: Bridge party
  • bubble_sort.pi: Bubble sort
  • carroll_symbolic_logic_60.pi: Lewis Carroll's symbolic logic problem #60 @@ -591,7 +600,7 @@

    Miscellaneous programs

  • josephus.pi: Josephus problem
  • lace.pi: Lace problem (ported from B-Prolog's entry to 2012 Prolog Programming Contest)
  • lcs_tabling.pi: Longest Common Subsequence using tabling (ported from B-Prolog) -
  • lunc_problem.pi: Lunch problem +
  • lunch_problem.pi: Lunch problem
  • mankell.pi: Generate all possible (right and wrong) spellings of 'Henning Mankell' and 'Kjellerstrand' (using a "grammar" etc)
  • matrix_slice.pi: Some list/array/matrix extraction functions. E.g. M.slice(1..3, 2..4) for extracting a sub square
  • mcsam.pi: MicroSAM program (ported from Sterling "The Art of Prolog") diff --git a/picat/knapsack_problem_mip.pi b/picat/knapsack_problem_mip.pi new file mode 100644 index 00000000..bb1a2059 --- /dev/null +++ b/picat/knapsack_problem_mip.pi @@ -0,0 +1,90 @@ +/* + + Knapsack problem (Rosetta Code) in Picat. + + From + http://rosettacode.org/wiki/Knapsack_problem/Unbounded + """ + A traveller gets diverted and has to make an unscheduled stop in what turns out to + be Shangri La. Opting to leave, he is allowed to take as much as he likes of the + following items, so long as it will fit in his knapsack, and he can carry it. He knows + that he can carry no more than 25 'weights' in total; and that the capacity of his + knapsack is 0.25 'cubic lengths'. + + Looking just above the bar codes on the items he finds their weights and volumes. + He digs out his recent copy of a financial paper and gets the value of each item. + + Item Explanation Value (each) weight Volume (each) + ------------------------------------------------------------ + panacea (vials of) Incredible healing properties 3000 0.3 0.025 + ichor (ampules of) Vampires blood 1800 0.2 0.015 + gold (bars) Shiney shiney 2500 2.0 0.002 + ---------------------------------------------------------------------------- + Knapsack For the carrying of - <=25 <=0.25 + + He can only take whole units of any item,, but there is much more of any item than + he could ever carry + + How many of each item does he take to maximise the value of items he is carrying + away with him? + + Note: + + 1. There are four solutions that maximise the value taken. Only one need be given. + """ + + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +go => + Value = [3000.0, 1800.0, 2500.0 ], + Weight = [ 0.3, 0.2, 2.0 ], + Volume = [ 0.025, 0.015, 0.002], + + N = Value.length, + + X = new_list(N), + X :: 0..1000, + + Z #= sum([X[I]*Value[I] : I in 1..N]), + + foreach(I in 1..N) + X[I] #>= 0 + end, + + knapsack(Volume, X, 0.25), + knapsack(Weight, X, 25.0), + + solve($[max(Z)], X), + + println(z=Z), + println(x=X), + Items = ["panacea","ichor","gold"], + foreach({Item,Num} in zip(Items,X), Num > 0) + printf("Take %d of %w\n", Num,Item) + end, + + println("\nTotal volume:"), + println(sum([X[I]*Volume[I] : I in 1..N])), + + println("Total weight:"), + println(sum([X[I]*Weight[I] : I in 1..N])), + + println("Total cost:"), + println(sum([X[I]*Value[I] : I in 1..N])), + + nl. + + +knapsack(W, Take, WTMax) => + sum([W[I]*Take[I] : I in 1..W.length]) #<= WTMax. \ No newline at end of file diff --git a/picat/least_square_optimeringslara_286.pi b/picat/least_square_optimeringslara_286.pi new file mode 100644 index 00000000..51fd8509 --- /dev/null +++ b/picat/least_square_optimeringslara_286.pi @@ -0,0 +1,47 @@ +/* + + Least square in Picat. + + Least square problem from the Swedish book + Optimeringslara, page 286f + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + +go => + + % temperature + T = [20.0, 30.0, 80.0,125.0,175.0,225.0,275.0,325.0,360.0,420.0,495.0,540.0,630.0,700.0], + % percentage gas + F = [0.0,5.8 ,14.7,31.6,43.2,58.3,78.4,89.4,96.4,99.1,99.5,99.9,100.0,100.0], + N = T.length, + + M = 4, + A = new_list(M+1), % equation (index 0..M) + A :: -100.0..100.0, + + Z :: -1000.0..1000.0, + Z #= sum([F[I] - (A[0+1] + A[1+1]*T[I] + A[2+1]*T[I]*T[I] + A[3+1]*T[I]*T[I]*T[I] + A[4+1]*T[I]*T[I]*T[I]*T[I]) : I in 1..N]), + + + A[0+1] + 20.0*A[1+1] + 20.0*20.0*A[2+1] + 20.0*20.0*20.0*A[3+1] + 20.0*20.0*20.0*20.0*A[4+1] #= 0.0, + A[0+1] + 700.0*A[1+1] + 700.0*700.0*A[2+1] + 700.0*700.0*700.0*A[3+1] + 700.0*700.0*700.0*700.0*A[4+1] #= 100.0, + foreach(I in 1..N) + (A[1+1] + 2.0*A[2+1]*T[I] + 3.0*A[3+1]*T[I]*T[I] + 4.0*A[4+1]*T[I]*T[I]*T[I]) #>= 0.0 + end, + + solve($[min(Z)],A), + + println(z=Z), + println(a=A), + + nl. + diff --git a/picat/markov_chains.pi b/picat/markov_chains.pi new file mode 100644 index 00000000..bda16958 --- /dev/null +++ b/picat/markov_chains.pi @@ -0,0 +1,122 @@ +/* + + Markov chains in Picat. + + From the (Swedish) book "Statistisk Dataanalys", page 299ff. + + Transition matrix of market shares for three products A, B, and C. + + From To + A B C + A 0.7 0.1 0,2 + B 0.2 0.6 0.2 + C 0.4 0.1 0.5 + + What is the stable state of the transitions? + + Also, this model solves the reversed problem: + given a stable state, generate a transition matrix. + + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + +go => + transitions(1,Transitions), + markov_chain(Transitions, X), + println(X), + nl. + +go2 => + transitions(2,Transitions), + markov_chain(Transitions, X), + println(X), + nl. + + +% reversed problem +go3 => + % This is the result from go/0 + X = [ 0.514285714285714, 0.2, 0.285714285714286], + markov_chain_rev(X, Transitions), + foreach(Row in Transitions) println(Row) end, + nl. + +% reversed problem 2 +go4 => + % This is the result from go2/0 + X = [0.367182763904145,0.07548189012462,0.102060886308352,0.15336027885601,0.034415433900105,0.080532266138785,0.13700725112718,0.023360147677889,0.01100743321783,0.015591648745084], + markov_chain_rev(X, Transitions), + + foreach(Row in Transitions) println(Row) end, + nl. + + + +% +% Given a transition matrix, yield stable state X +% +markov_chain(Transitions, X) => + + N = Transitions.length, + + X = new_list(N), + X :: 0.0..1.0, + + foreach(I in 1..N) + X[I] #= sum([Transitions[I,J]*X[J] : J in 1..N]) + end, + + sum(X) #= 1.0, + + solve(X). + + +transitions(1,Transitions) => + Transitions = + [[0.7,0.2,0.4], + [0.1,0.6,0.1], + [0.2,0.2,0.5]]. + +% Another problem from the book "Optimieringslara": +transitions(2,Transitions) => + Transitions = +[[0.01782,0.54024,0.46916,0.85467,0.30838,0.62996,0.42815,0.22189,0.37009,0.74681], + [0.04633,0.04427,0.14569,0.03971,0.50937,0.09334,0.01055,0.12776,0.36217,0.04528], + [0.12289,0.22650,0.12335,0.02186,0.10853,0.14409,0.05414,0.01577,0.01709,0.03761], + [0.28441,0.01805,0.09189,0.04428,0.00565,0.00538,0.16482,0.20185,0.15423,0.11391], + [0.00980,0.07623,0.03248,0.00753,0.01028,0.00424,0.12144,0.11167,0.05904,0.00009], + [0.15721,0.09004,0.07690,0.00448,0.03511,0.05270,0.00859,0.00057,0.00754,0.04810], + [0.29881,0.00202,0.00490,0.00909,0.00343,0.00296,0.17159,0.05433,0.00924,0.00041], + [0.02110,0.00004,0.02863,0.00255,0.01232,0.06338,0.02043,0.16868,0.00080,0.00129], + [0.01843,0.00002,0.02181,0.00453,0.00290,0.00011,0.00819,0.00118,0.00312,0.00163], + [0.02320,0.00259,0.00519,0.01130,0.00403,0.00384,0.01210,0.09630,0.01668,0.00487]]. + + +% +% reversed problem: Given stable state X, yield a transition matrix. +% +markov_chain_rev(X, Transitions) => + + N = length(X), + + Transitions = new_array(N,N), + Transitions.vars() :: 0.0..1.0, + + % Transitions[1,1] #<= 0.99, % to make it more interesting... + foreach(I in 1..N) + X[I] #= sum([Transitions[I,J]*X[J] : J in 1..N]) + end, + + sum(X) #= 1.0, + + + solve(Transitions). diff --git a/picat/markov_chains_taha.pi b/picat/markov_chains_taha.pi new file mode 100644 index 00000000..84ed0d68 --- /dev/null +++ b/picat/markov_chains_taha.pi @@ -0,0 +1,75 @@ +/* + + Markov chains in Picat. + + From Hamdy Taha "Operations Research" (8th edition), page 649ff. + Fertilizer example. + + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +go => + Cost = [100.0, 125.0, 160.0], + Mat = [[0.3, 0.6, 0.1], + [0.1, 0.6, 0.3], + [0.05, 0.4, 0.55]], + + % the transition matrix page 650 + % Mat = [[0.35, 0.6, 0.05], + % [0.3 , 0.6, 0.1 ], + % [0.25, 0.4, 0.35]], + + + N = Cost.length, + + %% Not used since it involves nonlinear * + % MeanFirstReturnTime = new_list(N), + % MeanFirstReturnTime :: 0.0..1.0, + + P = new_list(N), + P :: 0.0..1.0, + + TotCost #= sum([Cost[I]*P[I] : I in 1..N]), + + steady_state_prob(Mat, P), + % get_mean_first_return_time(P, MeanFirstReturnTime), + + % solve(P ++ MeanFirstReturnTime), + solve($[min(TotCost)], P), + + println(p=P), + println(totCost=TotCost), + % println(meanFirstReturnTime=MeanFirstReturnTime), + + nl. + + +% +% Calculates the steady state probablity of a transition matrix m +% +steady_state_prob(M, Prob) => + Len = M.length, + foreach(I in 1..Len) + Prob[I] #= sum([Prob[J]* M[J,I] : J in 1..Len]) + end, + sum(Prob) #= 1.0. + + +% +% calculate the mean first return time from a steady state probability array +% +% Gives error: nonlinear constraint,* +get_mean_first_return_time(Prob, Mfrt) => + foreach(I in 1..Prob.length) + Mfrt[I] * Prob[I] #= 1.0 + end. diff --git a/picat/mortgage.pi b/picat/mortgage.pi new file mode 100644 index 00000000..aed6dc6f --- /dev/null +++ b/picat/mortgage.pi @@ -0,0 +1,59 @@ +/* + + Mortgage in Picat. + + Marriot & Stuckey: Programming with Constraints, page 175f + (famous early example) + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +% b=203.128769950000162 +go => + mortgage(1000.0,10,10.0/100.0,150.0,B), + println(b=B), + nl. + +% p=921.685065855701964 +go2 => + mortgage(P,10,10.0/100.0,150.0,0), + println(p=P), + nl. + +% This give p=0.0,r=0.0,b=0.0 +% (book example: 0.3855*B + 6.1446*R) +go3 => + mortgage(P,10,10.0/100.0,R,B), + println([p=P,r=R,b=B]), + nl. + +go4 => + B #= 10.0, + mortgage(1234,10,10.0/100.0,R,B), + println([r=R,b=B]), + nl. + +% P: principal (amount borrowed) +% T: time periods +% I: interest rate +% R: repayment +% B: balance (final amount owing) +mortgage(P,T,I,R,B) ?=> + T = 0, + B #= P, + solve([P,T,I,R,B]). + +mortgage(P,T,I,R,B) => + T >= 0, + NT #= T -1, + NP #= P + P*I -R, + mortgage(NP,NT,I,R,B). diff --git a/picat/plan_mip.pi b/picat/plan_mip.pi new file mode 100644 index 00000000..cbf47a6a --- /dev/null +++ b/picat/plan_mip.pi @@ -0,0 +1,65 @@ +/* + + Planning model (MIP) in Picat. + + From GLPK:s example plan.mod + + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +go => + Bin1 :: 0.0..200.0, + Bin2 :: 0.0..2500.0, + Bin3 :: 400.0..800.0, + Bin4 :: 100.0..700.0, + Bin5 :: 0.0..1500.0, + Alum :: 0.0..1000.0, + Silicon :: 0.0..1000.0, + + Value #= 0.03 * Bin1 + 0.08 * Bin2 + 0.17 * Bin3 + 0.12 * Bin4 + 0.15 * Bin5 + 0.21 * Alum + 0.38 * Silicon, + + Alum #>= 0.0, + Silicon #>= 0.0, + Bin1 + Bin2 + Bin3 + Bin4 + Bin5 + Alum + Silicon #= 2000.0, + 0.15 * Bin1 + 0.04 * Bin2 + 0.02 * Bin3 + 0.04 * Bin4 + 0.02 * Bin5 + + 0.01 * Alum + 0.03 * Silicon #<= 60.0, + + 0.03 * Bin1 + 0.05 * Bin2 + 0.08 * Bin3 + 0.02 * Bin4 + 0.06 * Bin5 + + 0.01 * Alum #<= 100.0, + + 0.02 * Bin1 + 0.04 * Bin2 + 0.01 * Bin3 + 0.02 * Bin4 + 0.02 * Bin5 #<= 40.0, + + 0.02 * Bin1 + 0.03 * Bin2 + 0.01 * Bin5 #<= 30.0, + + 0.70 * Bin1 + 0.75 * Bin2 + 0.80 * Bin3 + 0.75 * Bin4 + 0.80 * Bin5 + 0.97 * Alum #>= 1500.0, + + 250.0 #<= 0.02 * Bin1 + 0.06 * Bin2 + 0.08 * Bin3 + 0.12 * Bin4 + 0.02 * Bin5 + 0.01 * Alum + 0.97 * Silicon, + + 0.02 * Bin1 + 0.06 * Bin2 + 0.08 * Bin3 + 0.12 * Bin4 + 0.02 * Bin5 + + 0.01 * Alum + 0.97 * Silicon #<= 300.0, + + + Vars = [Bin1,Bin2,Bin3,Bin4,Bin5,Alum,Silicon], + solve($[min(Value)], Vars), + + println(bin1=Bin1), + println(bin2=Bin2), + println(bin3=Bin3), + println(bin4=Bin4), + println(bin5=Bin5), + println(alum=Alum), + println(silicon=Silicon), + println(value=Value), + + nl. + diff --git a/picat/spreadsheet.pi b/picat/spreadsheet.pi new file mode 100644 index 00000000..d3f8df30 --- /dev/null +++ b/picat/spreadsheet.pi @@ -0,0 +1,41 @@ +/* + + Spreadsheet in Picat. + + From Krzysztof Apt "Principles of Constraint Programming" page 16ff. Spreadsheet. + Cf Winston "Artificial Intelligence", 3rd edition, page 235 + (not the same values though) + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +go => + + X = [B1,B4,B5,C4,C5,D4,D5,E7,E8], + X :: 0.0..1000.0, + + B1 #= 0.17, + B4 #= 3.5, + B5 #= 1.7, + C4 #= 1.5, + C5 #= 4.5, + D4 #= B4 * C4, + D5 #= B5 * C5, + E7 #= D4 + D5, + E8 #= E7 * (1.0 + B1), + + solve(X), + + println(X), + + nl. + diff --git a/picat/stigler.pi b/picat/stigler.pi new file mode 100644 index 00000000..40aec2e8 --- /dev/null +++ b/picat/stigler.pi @@ -0,0 +1,374 @@ +/* + + Original Stigler's 1939 diet problem in Picat. + + From GLPK:s example stigler.mod + """ + STIGLER, original Stigler's 1939 diet problem + + The Stigler Diet is an optimization problem named for George Stigler, + a 1982 Nobel Laureate in economics, who posed the following problem: + For a moderately active man weighing 154 pounds, how much of each of + 77 foods should be eaten on a daily basis so that the man's intake of + nine nutrients will be at least equal to the recommended dietary + allowances (RDSs) suggested by the National Research Council in 1943, + with the cost of the diet being minimal? + + The nutrient RDAs required to be met in Stigler's experiment were + calories, protein, calcium, iron, vitamin A, thiamine, riboflavin, + niacin, and ascorbic acid. The result was an annual budget allocated + to foods such as evaporated milk, cabbage, dried navy beans, and beef + liver at a cost of approximately $0.11 a day in 1939 U.S. dollars. + + While the name "Stigler Diet" was applied after the experiment by + outsiders, according to Stigler, "No one recommends these diets for + anyone, let alone everyone." The Stigler diet has been much ridiculed + for its lack of variety and palatability, however his methodology has + received praise and is considered to be some of the earliest work in + linear programming. + + The Stigler diet question is a linear programming problem. Lacking + any sophisticated method of solving such a problem, Stigler was + forced to utilize heuristic methods in order to find a solution. The + diet question originally asked in which quantities a 154 pound male + would have to consume 77 different foods in order to fulfill the + recommended intake of 9 different nutrients while keeping expense at + a minimum. Through "trial and error, mathematical insight and + agility," Stigler was able to eliminate 62 of the foods from the + original 77 (these foods were removed based because they lacked + nutrients in comparison to the remaining 15). From the reduced list, + Stigler calculated the required amounts of each of the remaining 15 + foods to arrive at a cost-minimizing solution to his question. + According to Stigler's calculations, the annual cost of his solution + was $39.93 in 1939 dollars. When corrected for inflation using the + consumer price index, the cost of the diet in 2005 dollars is + $561.43. The specific combination of foods and quantities is as + follows: + + Stigler's 1939 Diet + + Food Annual Quantities Annual Cost + ---------------- ----------------- ----------- + Wheat Flour 370 lb. $13.33 + Evaporated Milk 57 cans 3.84 + Cabbage 111 lb. 4.11 + Spinach 23 lb. 1.85 + Dried Navy Beans 285 lb. 16.80 + ---------------------------------------------- + Total Annual Cost $39.93 + + The 9 nutrients that Stigler's diet took into consideration and their + respective recommended daily amounts were: + + Table of nutrients considered in Stigler's diet + + Nutrient Daily Recommended Intake + ------------------------- ------------------------ + Calories 3,000 Calories + Protein 70 grams + Calcium .8 grams + Iron 12 milligrams + Vitamin A 5,000 IU + Thiamine (Vitamin B1) 1.8 milligrams + Riboflavin (Vitamin B2) 2.7 milligrams + Niacin 18 milligrams + Ascorbic Acid (Vitamin C) 75 milligrams + + Seven years after Stigler made his initial estimates, the development + of George Dantzig's Simplex algorithm made it possible to solve the + problem without relying on heuristic methods. The exact value was + determined to be $39.69 (using the original 1939 data). Dantzig's + algorithm describes a method of traversing the vertices of a polytope + of N+1 dimensions in order to find the optimal solution to a specific + situation. + + (From Wikipedia, the free encyclopedia.) + + Translated from GAMS by Andrew Makhorin . + + For the original GAMS model stigler1939.gms see [3]. + + References: + + 1. George J. Stigler, "The Cost of Subsistence," J. Farm Econ. 27, + 1945, pp. 303-14. + + 2. National Research Council, "Recommended Daily Allowances," Reprint + and Circular Series No. 115, January, 1943. + + 3. Erwin Kalvelagen, "Model building with GAMS," Chapter 2, "Building + linear programming models," pp. 128-34. + """ + + + Solution according GLPK: + """ + + Commodity Unit Quantity Cost + ------------------------- ---------- ---------- ---------- + Wheat Flour (Enriched) 10 lb. 29.95 $ 10.78 + Liver (Beef) 1 lb. 2.58 $ 0.69 + Cabbage 1 lb. 110.70 $ 4.10 + Spinach 1 lb. 22.58 $ 1.83 + Navy Beans, Dried 1 lb. 377.81 $ 22.29 + ----------------- + Total: $ 39.69 + """ + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +% import util. +import mip. + + +main => go. + + +go => + + stigler, + nl. + + +stigler => + data(Data), + NumCommodities = Data.length, + + % days in a year + Days = 365.25, + + % nutrients + NumNutrients = 9, + + allowance(Allowance), + + % decision variables + + % dollars of food to be purchased daily + X = new_list(NumCommodities), + X :: 0.0..1.0, + + XCost = new_list(NumCommodities), + XCost :: 0.0..100.0, + Quant = new_list(NumCommodities), + Quant :: 0.0..1000.0, + + % total food bill + TotalCost #= Days * sum(X), + Cost #= sum(X), % cost per day + + % nutrient balance + foreach(N in 3..NumNutrients+2) + sum([Data[C,N] * X[C] : C in 1..NumCommodities]) #>= Allowance[N-2] + end, + + foreach(C in 1..NumCommodities) + XCost[C] #= Days * X[C], + % Quant[C] #= 100.0*Days*X[C] / Data[C,1] % / is integer divistion + Quant[C]*Data[C,1] #= 100.0*Days*X[C] + end, + + Vars = X ++ XCost ++ Quant, + solve($[min(Cost)], Vars), + + % println(x=X), + println(cost=Cost), + % println(xcost=XCost), + % println(quant=Quant), + println(totalCost=TotalCost), + + commodities(Commodities), + foreach({Comm,Q,C} in zip(Commodities,Quant,XCost), Q > 0) + printf("%-35w: %10.2f $ %5.2f\n", Comm,Q,C) + end, + + nl. + +allowance(Allowance) => + Allowance = [3.0, 70.0, 0.8, 12.0, 5.0, 1.8, 2.7, 18.0, 75.0]. + + +commodities(Commodities) => + Commodities = +[ +"Wheat Flour (Enriched) 10 lb.", +"Macaroni 1 lb.", +"Wheat Cereal (Enriched) 28 oz.", +"Corn Flakes 8 oz.", +"Corn Meal 1 lb.", +"Hominy Grits 24 oz.", +"Rice 1 lb.", +"Rolled Oats 1 lb.", +"White Bread (Enriched) 1 lb.", +"Whole Wheat Bread 1 lb.", +"Rye Bread 1 lb.", +"Pound Cake 1 lb.", +"Soda Crackers 1 lb.", +"Milk 1 qt.", +"Evaporated Milk (can) 14.5 oz.", +"Butter 1 lb.", +"Oleomargarine 1 lb.", +"Eggs 1 doz.", +"Cheese (Cheddar) 1 lb.", +"Cream 1/2 pt.", +"Peanut Butter 1 lb.", +"Mayonnaise 1/2 pt.", +"Crisco 1 lb.", +"Lard 1 lb.", +"Sirloin Steak 1 lb.", +"Round Steak 1 lb.", +"Rib Roast 1 lb.", +"Chuck Roast 1 lb.", +"Plate 1 lb.", +"Liver (Beef) 1 lb.", +"Leg of Lamb 1 lb.", +"Lamb Chops (Rib) 1 lb.", +"Pork Chops 1 lb.", +"Pork Loin Roast 1 lb.", +"Bacon 1 lb.", +"Ham - smoked 1 lb.", +"Salt Pork 1 lb.", +"Roasting Chicken 1 lb.", +"Veal Cutlets 1 lb.", +"Salmon, Pink (can) 16 oz.", +"Apples 1 lb.", +"Bananas 1 lb.", +"Lemons 1 doz.", +"Oranges 1 doz.", +"Green Beans 1 lb.", +"Cabbage 1 lb.", +"Carrots 1 bunch", +"Celery 1 stalk", +"Lettuce 1 head", +"Onions 1 lb.", +"Potatoes 15 lb.", +"Spinach 1 lb.", +"Sweet Potatoes 1 lb.", +"Peaches (can) No. 2 1/2", +"Pears (can) No. 2 1/2", +"Pineapple (can) No. 2 1/2", +"Asparagus (can) No. 2", +"Grean Beans (can) No. 2", +"Pork and Beans (can) 16 oz.", +"Corn (can) No. 2", +"Peas (can) No. 2", +"Tomatoes (can) No. 2", +"Tomato Soup (can) 10 1/2 oz.", +"Peaches, Dried 1 lb.", +"Prunes, Dried 1 lb.", +"Raisins, Dried 15 oz.", +"Peas, Dried 1 lb.", +"Lima Beans, Dried 1 lb.", +"Navy Beans, Dried 1 lb.", +"Coffee 1 lb.", +"Tea 1/4 lb.", +"Cocoa 8 oz.", +"Chocolate 8 oz.", +"Sugar 10 lb.", +"Corn Sirup 24 oz.", +"Molasses 18 oz.", +"Strawberry Preserve 1 lb."]. + + + +% set N := +% calories % Calories, unit = 1000 +% protein % Protein, unit = grams +% calcium % Calcium, unit = grams +% iron % Iron, unit = milligrams +% vitaminA % Vitamin A, unit = 1000 International Units +% thiamine % Thiamine, Vit. B1, unit = milligrams +% riboflavin % Riboflavin, Vit. B2, unit = milligrams +% niacin % Niacin (Nicotinic Acid), unit = milligrams +% ascorbicAcid % Ascorbic Acid, Vit. C, unit = milligrams +% ; + +% price and weight is the two first columns +data(Data) => + % 1..num_nutrients+2, +Data = +[ + [36.0, 12600.0, 44.7, 1411.0, 2.0, 365.0, 0.0, 55.4, 33.3, 441.0, 0.0], + [14.1, 3217.0, 11.6, 418.0, 0.7, 54.0, 0.0, 3.2, 1.9, 68.0, 0.0], + [24.2, 3280.0, 11.8, 377.0, 14.4, 175.0, 0.0, 14.4, 8.8, 114.0, 0.0], + [ 7.1, 3194.0, 11.4, 252.0, 0.1, 56.0, 0.0, 13.5, 2.3, 68.0, 0.0], + [ 4.6, 9861.0, 36.0, 897.0, 1.7, 99.0, 30.9, 17.4, 7.9, 106.0, 0.0], + [ 8.5, 8005.0, 28.6, 680.0, 0.8, 80.0, 0.0, 10.6, 1.6, 110.0, 0.0], + [ 7.5, 6048.0, 21.2, 460.0, 0.6, 41.0, 0.0, 2.0, 4.8, 60.0, 0.0], + [ 7.1, 6389.0, 25.3, 907.0, 5.1, 341.0, 0.0, 37.1, 8.9, 64.0, 0.0], + [ 7.9, 5742.0, 15.6, 488.0, 2.5, 115.0, 0.0, 13.8, 8.5, 126.0, 0.0], + [ 9.1, 4985.0, 12.2, 484.0, 2.7, 125.0, 0.0, 13.9, 6.4, 160.0, 0.0], + [ 9.2, 4930.0, 12.4, 439.0, 1.1, 82.0, 0.0, 9.9, 3.0, 66.0, 0.0], + [24.8, 1829.0, 8.0, 130.0, 0.4, 31.0, 18.9, 2.8, 3.0, 17.0, 0.0], + [15.1, 3004.0, 12.5, 288.0, 0.5, 50.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [11.0, 8867.0, 6.1, 310.0, 10.5, 18.0, 16.8, 4.0, 16.0, 7.0, 177.0], + [ 6.7, 6035.0, 8.4, 422.0, 15.1, 9.0, 26.0, 3.0, 23.5, 11.0, 60.0], + [20.8, 1473.0, 10.8, 9.0, 0.2, 3.0, 44.2, 0.0, 0.2, 2.0, 0.0], + [16.1, 2817.0, 20.6, 17.0, 0.6, 6.0, 55.8, 0.2, 0.0, 0.0, 0.0], + [32.6, 1857.0, 2.9, 238.0, 1.0, 52.0, 18.6, 2.8, 6.5, 1.0, 0.0], + [24.2, 1874.0, 7.4, 448.0, 16.4, 19.0, 28.1, 0.8, 10.3, 4.0, 0.0], + [14.1, 1689.0, 3.5, 49.0, 1.7, 3.0, 16.9, 0.6, 2.5, 0.0, 17.0], + [17.9, 2534.0, 15.7, 661.0, 1.0, 48.0, 0.0, 9.6, 8.1, 471.0, 0.0], + [16.7, 1198.0, 8.6, 18.0, 0.2, 8.0, 2.7, 0.4, 0.5, 0.0, 0.0], + [20.3, 2234.0, 20.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [ 9.8, 4628.0, 41.7, 0.0, 0.0, 0.0, 0.2, 0.0, 0.5, 5.0, 0.0], + [39.6, 1145.0, 2.9, 166.0, 0.1, 34.0, 0.2, 2.1, 2.9, 69.0, 0.0], + [36.4, 1246.0, 2.2, 214.0, 0.1, 32.0, 0.4, 2.5, 2.4, 87.0, 0.0], + [29.2, 1553.0, 3.4, 213.0, 0.1, 33.0, 0.0, 0.0, 2.0, 0.0, 0.0], + [22.6, 2007.0, 3.6, 309.0, 0.2, 46.0, 0.4, 1.0, 4.0, 120.0, 0.0], + [14.6, 3107.0, 8.5, 404.0, 0.2, 62.0, 0.0, 0.9, 0.0, 0.0, 0.0], + [26.8, 1692.0, 2.2, 333.0, 0.2, 139.0, 169.2, 6.4, 50.8, 316.0, 525.0], + [27.6, 1643.0, 3.1, 245.0, 0.1, 20.0, 0.0, 2.8, 3.0, 86.0, 0.0], + [36.6, 1239.0, 3.3, 140.0, 0.1, 15.0, 0.0, 1.7, 2.7, 54.0, 0.0], + [30.7, 1477.0, 3.5, 196.0, 0.2, 80.0, 0.0, 17.4, 2.7, 60.0, 0.0], + [24.2, 1874.0, 4.4, 249.0, 0.3, 37.0, 0.0, 18.2, 3.6, 79.0, 0.0], + [25.6, 1772.0, 10.4, 152.0, 0.2, 23.0, 0.0, 1.8, 1.8, 71.0, 0.0], + [27.4, 1655.0, 6.7, 212.0, 0.2, 31.0, 0.0, 9.9, 3.3, 50.0, 0.0], + [16.0, 2835.0, 18.8, 164.0, 0.1, 26.0, 0.0, 1.4, 1.8, 0.0, 0.0], + [30.3, 1497.0, 1.8, 184.0, 0.1, 30.0, 0.1, 0.9, 1.8, 68.0, 46.0], + [42.3, 1072.0, 1.7, 156.0, 0.1, 24.0, 0.0, 1.4, 2.4, 57.0, 0.0], + [13.0, 3489.0, 5.8, 705.0, 6.8, 45.0, 3.5, 1.0, 4.9, 209.0, 0.0], + [ 4.4, 9072.0, 5.8, 27.0, 0.5, 36.0, 7.3, 3.6, 2.7, 5.0, 544.0], + [ 6.1, 4982.0, 4.9, 60.0, 0.4, 30.0, 17.4, 2.5, 3.5, 28.0, 498.0], + [26.0, 2380.0, 1.0, 21.0, 0.5, 14.0, 0.0, 0.5, 0.0, 4.0, 952.0], + [30.9, 4439.0, 2.2, 40.0, 1.1, 18.0, 11.1, 3.6, 1.3, 10.0, 1993.0], + [ 7.1, 5750.0, 2.4, 138.0, 3.7, 80.0, 69.0, 4.3, 5.8, 37.0, 862.0], + [ 3.7, 8949.0, 2.6, 125.0, 4.0, 36.0, 7.2, 9.0, 4.5, 26.0, 5369.0], + [ 4.7, 6080.0, 2.7, 73.0, 2.8, 43.0, 188.5, 6.1, 4.3, 89.0, 608.0], + [ 7.3, 3915.0, 0.9, 51.0, 3.0, 23.0, 0.9, 1.4, 1.4, 9.0, 313.0], + [ 8.2, 2247.0, 0.4, 27.0, 1.1, 22.0, 112.4, 1.8, 3.4, 11.0, 449.0], + [ 3.6, 11844.0, 5.8, 166.0, 3.8, 59.0, 16.6, 4.7, 5.9, 21.0, 1184.0], + [34.0, 16810.0, 14.3, 336.0, 1.8, 118.0, 6.7, 29.4, 7.1, 198.0, 2522.0], + [ 8.1, 4592.0, 1.1, 106.0, 0.0, 138.0, 918.4, 5.7, 13.8, 33.0, 2755.0], + [ 5.1, 7649.0, 9.6, 138.0, 2.7, 54.0, 290.7, 8.4, 5.4, 83.0, 1912.0], + [16.8, 4894.0, 3.7, 20.0, 0.4, 10.0, 21.5, 0.5, 1.0, 31.0, 196.0], + [20.4, 4030.0, 3.0, 8.0, 0.3, 8.0, 0.8, 0.8, 0.8, 5.0, 81.0], + [21.3, 3993.0, 2.4, 16.0, 0.4, 8.0, 2.0, 2.8, 0.8, 7.0, 399.0], + [27.7, 1945.0, 0.4, 33.0, 0.3, 12.0, 16.3, 1.4, 2.1, 17.0, 272.0], + [10.0, 5386.0, 1.0, 54.0, 2.0, 65.0, 53.9, 1.6, 4.3, 32.0, 431.0], + [ 7.1, 6389.0, 7.5, 364.0, 4.0, 134.0, 3.5, 8.3, 7.7, 56.0, 0.0], + [10.4, 5452.0, 5.2, 136.0, 0.2, 16.0, 12.0, 1.6, 2.7, 42.0, 218.0], + [13.8, 4109.0, 2.3, 136.0, 0.6, 45.0, 34.9, 4.9, 2.5, 37.0, 370.0], + [ 8.6, 6263.0, 1.3, 63.0, 0.7, 38.0, 53.2, 3.4, 2.5, 36.0, 1253.0], + [ 7.6, 3917.0, 1.6, 71.0, 0.6, 43.0, 57.9, 3.5, 2.4, 67.0, 862.0], + [15.7, 2889.0, 8.5, 87.0, 1.7, 173.0, 86.8, 1.2, 4.3, 55.0, 57.0], + [ 9.0, 4284.0, 12.8, 99.0, 2.5, 154.0, 85.7, 3.9, 4.3, 65.0, 257.0], + [ 9.4, 4524.0, 13.5, 104.0, 2.5, 136.0, 4.5, 6.3, 1.4, 24.0, 136.0], + [ 7.9, 5742.0, 20.0, 1367.0, 4.2, 345.0, 2.9, 28.7, 18.4, 162.0, 0.0], + [ 8.9, 5097.0, 17.4, 1055.0, 3.7, 459.0, 5.1, 26.9, 38.2, 93.0, 0.0], + [ 5.9, 7688.0, 26.9, 1691.0, 11.4, 792.0, 0.0, 38.4, 24.6, 217.0, 0.0], + [22.4, 2025.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 5.1, 50.0, 0.0], + [17.4, 652.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.3, 42.0, 0.0], + [ 8.6, 2637.0, 8.7, 237.0, 3.0, 72.0, 0.0, 2.0, 11.9, 40.0, 0.0], + [16.2, 1400.0, 8.0, 77.0, 1.3, 39.0, 0.0, 0.9, 3.4, 14.0, 0.0], + [51.7, 8773.0, 34.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [13.7, 4996.0, 14.7, 0.0, 0.5, 74.0, 0.0, 0.0, 0.0, 5.0, 0.0], + [13.6, 3752.0, 9.0, 0.0, 10.3, 244.0, 0.0, 1.9, 7.5, 146.0, 0.0], + [20.5, 2213.0, 6.4, 11.0, 0.4, 7.0, 0.2, 0.2, 0.4, 3.0, 0.0] +]. + + diff --git a/picat/transp.pi b/picat/transp.pi new file mode 100644 index 00000000..40c7d0f5 --- /dev/null +++ b/picat/transp.pi @@ -0,0 +1,92 @@ +/* + + Transportation problem in Picat. + + From GLPK:s example transp.mod + + """ + A TRANSPORTATION PROBLEM + + This problem finds a least cost shipping schedule that meets + requirements at markets and supplies at factories. + + References: + Dantzig G B, "Linear Programming and Extensions." + Princeton University Press, Princeton, New Jersey, 1963, + Chapter 3-3. + """ + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +go => + + % canning plants + % set I := Seattle San-Diego; + NumPlants = 2, + + % markets + % set J := New-York Chicago Topeka; + NumMarkets = 3, + + + % capacity of plant i in cases + A = [350.0, 600.0], + + % demand at market j in cases + B = [325.0, 300.0, 275.0], + + % distance in thousands of miles + D = [ + [2.5, 1.7, 1.8], + [2.5, 1.8, 1.4] + ], + + % freight in dollars per case per thousand miles + F = 90.0, + + + % transport cost in thousands of dollars per case + % param c{i in I, j in J} := f * d[i,j] / 1000; + % array[I,J] of float: c = array2d(I,J, [f*d[i,j] / 1000.0 | i in I, j in J]); + C = [ [F*D[I,J] : J in 1..NumMarkets] : I in 1..NumPlants], + + + X = new_array(NumPlants,NumMarkets), + X.vars() :: 0.0..1000.0, + + Cost * 1000.0 #= sum([C[I,J]*X[I,J] : I in 1..NumPlants, J in 1..NumMarkets]), + + % observe supply limit at plant i + foreach(I in 1..NumPlants) + sum([X[I,J] : J in 1..NumMarkets]) #<= A[I] + end, + + + % satisfy demand at market j + foreach(J in 1..NumMarkets) + sum([X[I,J] : I in 1..NumPlants]) #>= B[J] + end, + + solve($[min(Cost)],X), + + printf("cost: %0.3f\n", Cost), + println("X:"), + foreach(Row in X) + foreach(R in Row) + printf("%8.2f ",R) + end, + nl + end, + + nl. + diff --git a/picat/transp1.pi b/picat/transp1.pi new file mode 100644 index 00000000..aed39656 --- /dev/null +++ b/picat/transp1.pi @@ -0,0 +1,136 @@ +/* + + Transportation problem in Picat. + + This is in princple the same model as in + Pascal Van Hentenryck "The OPL Optimization Programming Language", + page 116. + + Data from OPL example transp1.dat + + + This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com + See also my Picat page: http://www.hakank.org/picat/ + +*/ + + +import mip. + + +main => go. + + +go => + capacity(Capacity), + supply(Supply), + demand(Demand), + cost(Cost), + + NumProducts = Supply.length, + NumCities = Supply[1].length, + + % decision variables + Trans = new_array(NumProducts,NumCities,NumCities), + % Trans.vars() :: 0.0..1000000.0, + + % Z :: 0.0..10000000.0, + Z #= sum([ Cost[P,O,D] * Trans[P,O,D] : P in 1..NumProducts , O in 1..NumCities, D in 1..NumCities]), + + + foreach(P in 1..NumProducts , O in 1..NumCities) + sum([Trans[P,O,D] : D in 1..NumCities]) #= Supply[P,O] + end, + + foreach(P in 1..NumProducts , D in 1..NumCities) + sum([Trans[P,O,D] : O in 1..NumCities]) #= Demand[P,D] + end, + + foreach(O in 1..NumCities, D in 1..NumCities) + sum([Trans[P,O,D] : P in 1..NumProducts]) #<= Capacity + end, + + + solve($[min(Z)],Trans), + + println(z=Z), + println("Trans:"), + foreach(T in Trans) + foreach(R in T) + foreach(S in R) + printf("%6.1f ", S) + end, + nl + end + end, + + + nl. + + +capacity(Capacity) => Capacity = 625.0. + +% 1..NumProducts,1..NumCities +supply(Supply) => + Supply = + [ + [400.0, 700.0, 800.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [800.0, 1600.0, 1800.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [200.0, 300.0, 300.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + ]. + +% 1..NumProducts,1..NumCities +demand(Demand) => + Demand = + [ + [0, 0, 0,300,300,100, 75,650,225,250], % bands + [0, 0, 0,500,750,400,250,950,850,500], % coils + [0, 0, 0,100,100, 0, 50,200,100,250] % plate + ]. + +% 1..NumProducts, 1..NumCities,1..NumCities +cost(Cost) => + Cost = + [ + % bands + [ + [0, 0, 0, 30, 10, 8, 10, 11, 71, 6], + [0, 0, 0, 22, 7, 10, 7, 21, 82, 13], + [0, 0, 0, 19, 11, 12, 10, 25, 83, 15], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ] + , + % coils + [ + [0, 0, 0, 39, 14, 11, 14, 16, 82, 8], + [0, 0, 0, 27, 9, 12, 9, 26, 95, 17], + [0, 0, 0, 24, 14, 17, 13, 28, 99, 20], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ] + , + % plate + [ + [0, 0, 0, 41, 15, 12, 16, 17, 86, 8], + [0, 0, 0, 29, 9, 13, 9, 28, 99, 18], + [0, 0, 0, 26, 14, 17, 13, 31, 104, 20], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ] + ].