MODEL:
 ! Product cycling Problem. 
Given:
  a set of N products, each with a:
   demand rate D, constant over time, and
   production rate P, when the product is being produced,
   holding cost rate  H per unit in inventory,
 and for each pair of products, (i, j):
   T(i,j) = changeover time from product i to product j,
   S(i,j) = changeover cost from product i to product j,
 Find:
   a sequence of the products that is in fact a cycle that repeats,
   a cycle length CYCTIME, 
  so as to
   minimize the cost of changeover + inventory, subject
   to cycle being long enough to produce every product
   exactly once.

 ! Keywords: TSP, traveling sales person, routing, tour,
   lotsizing, cyclic production, Economic Manufacturing Quantity,
   EMQ;
SETS:
  PROD: D, P, H,
        U; ! U( I) = sequence no. of PROD;
  LINK( PROD, PROD):
       T,  ! Changeover time matrix;
       S,  ! Changeover cost matrix;
    X, Y;  ! X( I, J) = 1 if link I, J is in tour;
 ENDSETS
DATA: PROD = Vanilla Choco Strawb BPecan Neopol; D = 2 1.5 1.3 1.4 1.1; !Demand rate; P = 21 20 22 19 16; !Production rate; H = 1 1.5 1.8 1.2 1.9; !Holding cost rate; !Changeover times ...; T = 0 .2 .2 .4 .5 ! from Van; .8 0 .4 .3 .2 ! from Choco; .6 .1 0 .4 .3 ! from Strawb; .4 .1 .1 0 .2 ! from BPecan; .2 .4 .2 .3 0; ! From Neopol; S = 0 4 5 5 6 !Changeover costs; 9 0 7 6 5 7 4 0 5 6 6 4 5 0 7 5 4 9 8 0; ENDDATA !-----------------------------------------------------------; ! Variables: X(I,J) = 1 if product J follows product I in the cycle, CYCTIME = length of the cycle, U(J) = sequence number of product J in the cycle, arbitrarily, U(1) = 0, Y(I,J) = 1 if the sequence I -> J is on the path from 1 to N (the highest numbered product). ; ! Note, the demand for product I during a cycle is D(I)*CYCTIME. Thus, the production time that must be devoted to product I during the cycle is CYCTIME*( D(I)/P(I)). We assume the inventory of each product hits 0 just before its production starts. The maximum inventory for product I is when production ends for the product. Thus, the maximum inventory is CYCTIME*( D(I)/P(I))*(P(I) - D(I)). The average inventory is half the max. Minimize average cost per unit time of changeover cost + inventory cost.; MIN = CYCOST/CYCTIME + INVCOST; CYCOST = @SUM( LINK(I,J): S(I,J)*X(I,J)); INVCOST = @SUM( PROD( I): H(I)*D(I)*CYCTIME*(1-D(I)/P(I))/2); ! The production cycle must be long enough so that the setup times + the production times <= CYCTIME; @SUM( LINK(I,J): T(I,J)* X(I,J)) + CYCTIME *@SUM( PROD( I): D(I)/P(I)) <= CYCTIME ; ! Warning: May take long to solve cases with N > 12. Set N = number of products/flavors; N = @SIZE( PROD); ! The following constraints force each product to appear exactly once in the sequence; @FOR( PROD( K): ! It must be entered; @SUM( PROD( I)| I #NE# K: X( I, K)) = 1; ! It must be departed; @SUM( PROD( J)| J #NE# K: X( K, J)) = 1; X( K, K) = 0; ! A weak form of the subtour breaking constraints; ! Not very powerful for large(N>10) problems; @FOR( PROD( J)| J #GT# 1 #AND# J #NE# K: U( J) >= U( K) + X ( K, J) - ( N - 2) * ( 1 - X( K, J)) + ( N - 3) * X( J, K); ); ); ! Make the X's 0/1; @FOR( LINK: @BIN( X);); ! Some more cuts. They may make it solve faster. Can be removed to make model smaller; ! For the first and last stop we know...; @FOR( PROD( K)| K #GT# 1: U( K) <= N - 1 - ( N - 2) * X( 1, K); U( K) >= 1 + ( N - 2) * X( K, 1); ); ! We know the sum of the stop numbers; @SUM( PROD( I): U( I)) = ( N-1)*N/2; ! Add some 'multi-commodity cuts'. One could have a commodity for every PROD. Here we do it only for PROD N; ! There has to be a path from 1 to N. Y(I,J) = 1 if (I,J) is on the path from 1 to N. It helps if 1 and N are distant; ! The PROD N's commodity must leave 1; @SUM( PROD( J)| J #NE# 1: Y( 1, J)) = 1; ! It must get to N; @SUM( PROD( I)| I #NE# N: Y( I, N)) = 1; ! If it goes into K, it must depart K; @FOR( PROD( K)| 1 #NE# K #AND# K #NE# N: @SUM( PROD( I)|I #NE# N: Y( I, K)) = @SUM( PROD( J)| J #NE# 1: Y( K, J)); ); ! If (I,J) is on path from 1 to N, it must also be on the full tour; @FOR( LINK( I, J): Y( I, J) <= X( I, J);); END