/*****************************************************************/ /* */ /* Program to Solve the */ /* Multi-level Network Optimization Problem */ /* */ /* Lagrangean Relaxation Routine */ /* */ /* Author: Frederico R. B. da Cruz */ /* Department of Industrial Engineering and */ /* Operations Research */ /* University of Massachusetts at Amherst */ /* USA */ /* */ /* E-mail: fcruz@sml1.ecs.umass.edu */ /* */ /* Version: 3.00 */ /* */ /* Date: May/98 */ /* */ /*****************************************************************/ #include #include #include "sl.c" #include "sm.c" void Lagr_Relax (int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *DemaNodLev, ListNodPtrType *RootNodLev, MlType *Ml, slType *sl, PredLevPtrType **GoodPredLev, PredLevPtrType **CurrPredLev, float *BestUpper, float *UpperBound, float *LowerBound); void Bounds (int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *DemaNodLev, ListNodPtrType *RootNodLev, MlType *Ml, slType *sl, PredLevPtrType **GoodPredLev, PredLevPtrType **CurrPredLev, float *UpperBound, float *LowerBound); float ComputeSqrNormGamma(int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *RootNodLev, MlType *Ml, slType *sl); void First_Ck (int SizeOfProblem, int *Base, int *Displacement, int *k, float *Ck); void Next_Ck (int *Base, int *Displacement, int *k, float *Ck); void UpdateMultipliers (int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *RootNodLev, float tk); /*****************************************************************/ void Lagr_Relax (int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *DemaNodLev, ListNodPtrType *RootNodLev, MlType *Ml, slType *sl, PredLevPtrType **GoodPredLev, PredLevPtrType **CurrPredLev, float *BestUpper, float *UpperBound, float *LowerBound) { int i,k,l; static int SizeOfProblem = 0; static int FirstNode = YES; int SizeLevel; int Base, Displacement; float Solution_L, Solution_M, Gap, SqrNormOfGamma, tk, Ck; PredLevPtrType *AuxPredLev; ListArcPtrType ListArcPtr; ListNodPtrType ListNodPtr; if ( SizeOfProblem == 0 ) { fprintf(OUTFILE, "\nProblem Size =" ); for ( i = 0; i < n_Nodes; i++ ) { ListArcPtr = AdjaNodInd[i]; while ( ListArcPtr != NULL ) { SizeOfProblem++; ListArcPtr = ListArcPtr->Next; } } fprintf(OUTFILE, " %d * %d", SizeOfProblem, n_Levels ); SizeOfProblem *= n_Levels; for ( l = 0; l < n_Levels; l++ ) { SizeLevel = 0; ListNodPtr = RootNodLev[l]; while ( ListNodPtr != NULL) { SizeLevel++; ListNodPtr = ListNodPtr->Next; } if ( SizeLevel == 1 ) SizeLevel = 0; SizeOfProblem += SizeLevel; fprintf(OUTFILE, " + %d", SizeLevel); } fprintf(OUTFILE, " = %d\n", SizeOfProblem ); } *LowerBound = -INFINITY; *UpperBound = INFINITY; First_Ck(SizeOfProblem, &Base, &Displacement, &k, &Ck); do { /* compute lower bound */ Solution_L = Solve_L(n_Nodes, n_Levels, AdjaNodInd, DemaNodLev, RootNodLev, *CurrPredLev, Ml, sl); if ( Solution_L > *LowerBound ) *LowerBound = Solution_L; /* compute upper bound */ if ( Solution_L < INFINITY ) Solution_M = Solve_M(n_Nodes, n_Levels, AdjaNodInd, RootNodLev); else Solution_M = INFINITY; if ( Solution_M < *UpperBound ) { *UpperBound = Solution_M; AuxPredLev = *GoodPredLev; *GoodPredLev = *CurrPredLev; *CurrPredLev = AuxPredLev; } /* compute gap */ /* if ( *LowerBound > 0 ) Gap = *UpperBound / *LowerBound - 1; else if ( *LowerBound < 0 ) Gap = -( *UpperBound / *LowerBound - 1); else Gap = INFINITY; */ Gap = 1 - ( *LowerBound / *UpperBound ); /* compute square of subgradient Euclidean norm */ SqrNormOfGamma = ComputeSqrNormGamma(n_Nodes, n_Levels, AdjaNodInd, RootNodLev, Ml, sl); /* apply subgradient optimization algorithm */ if ( ( Gap > EPSILON ) && ( SqrNormOfGamma > ZERO ) ) { /* compute step */ tk = Ck * ( *UpperBound-Solution_L ) / SqrNormOfGamma; /* fprintf(OUTFILE, "\n k -> C_k * ( UB - Sol_k) / Sqr(Gamma_k) = Step\n"); fprintf(OUTFILE, "%3d -> %f * (%6.0f - %8.1f) / %10.2f = %f\n", k, Ck, *UpperBound, Solution_L, SqrNormOfGamma, tk); */ /* update Lagrangean multipliers */ UpdateMultipliers(n_Nodes, n_Levels, AdjaNodInd, RootNodLev, tk); /* update Ck */ Next_Ck(&Base, &Displacement, &k, &Ck); } } while ( ( *LowerBound < *BestUpper ) && ( Gap > EPSILON ) && ( SqrNormOfGamma > ZERO ) && #if _STEP_ == _HO_SE_ ( k <= MAXITERATIONS ) ); #elif _STEP_ == _HE_WO_CR_ ( tk > ZERO ) ); #endif /* write partial results */ if ( FirstNode == YES ) { FirstNode = NO; /* fprintf(OUTFILE, "\nFirst Node:"); */ fprintf(OUTFILE, "\nIterations %d", k); } } void Bounds (int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *DemaNodLev, ListNodPtrType *RootNodLev, MlType *Ml, slType *sl, PredLevPtrType **GoodPredLev, PredLevPtrType **CurrPredLev, float *UpperBound, float *LowerBound) { float Solution_L, Solution_M; PredLevPtrType *AuxPredLev; /* initialize bounds */ *LowerBound = -INFINITY; *UpperBound = INFINITY; /* compute lower bound */ Solution_L = Solve_L(n_Nodes, n_Levels, AdjaNodInd, DemaNodLev, RootNodLev, *CurrPredLev, Ml, sl); if ( Solution_L > *LowerBound ) *LowerBound = Solution_L; /* compute upper bound */ if ( Solution_L < INFINITY ) Solution_M = Solve_M(n_Nodes, n_Levels, AdjaNodInd, RootNodLev); else Solution_M = INFINITY; if ( Solution_M < *UpperBound ) { *UpperBound = Solution_M; AuxPredLev = *GoodPredLev; *GoodPredLev = *CurrPredLev; *CurrPredLev = AuxPredLev; } } float ComputeSqrNormGamma(int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *RootNodLev, MlType *Ml, slType *sl) { int i, l; float SqrSum; ListArcPtrType ListArcPtr; ListNodPtrType ListNodPtr; SqrSum = 0; for ( l =0; l < n_Levels; l++ ) { ListNodPtr = RootNodLev[l]; while ( ListNodPtr != NULL) { if ( ListNodPtr->NPtr->Status == SET_FREE ) { ListNodPtr->NPtr->Gamma_vi = ListNodPtr->NPtr->TransFlow - sl[l] * ListNodPtr->NPtr->zi; SqrSum += ListNodPtr->NPtr->Gamma_vi * ListNodPtr->NPtr->Gamma_vi; } ListNodPtr = ListNodPtr->Next; } for ( i = 0; i < n_Nodes; i++ ) { ListArcPtr = AdjaNodInd[i]; while ( ListArcPtr != NULL ) { if ( ListArcPtr->APtr[l]->Status == SET_FREE ) { ListArcPtr->APtr[l]->Gamma_wij = ListArcPtr->APtr[l]->xij - Ml[l] * ListArcPtr->APtr[l]->yij; SqrSum += ListArcPtr->APtr[l]->Gamma_wij * ListArcPtr->APtr[l]->Gamma_wij; } ListArcPtr = ListArcPtr->Next; } } } /* fprintf(OUTFILE,"\nGammas\n" ); fprintf(OUTFILE,"\n( i) -> ( TransFlow - sl * zi) = Gamma_vi\n"); for ( l =0; l < n_Levels; l++ ) { ListNodPtr = RootNodLev[l]; while ( ListNodPtr != NULL) { if ( ListNodPtr->NPtr->Status == SET_FREE ) fprintf(OUTFILE, "(%2d) -> ( %f - %f * %d) = %10.6f\n", ListNodPtr->NPtr->i+1, ListNodPtr->NPtr->TransFlow, sl[l], ListNodPtr->NPtr->zi, ListNodPtr->NPtr->Gamma_vi); ListNodPtr = ListNodPtr->Next; } } fprintf(OUTFILE, "\n( l, i, j) -> ( xlij - Ml * ylij) = Gamma_wlij\n" ); for (i=0; i < n_Nodes; i++ ) { for ( l =0; l < n_Levels; l++ ) { ListArcPtr = AdjaNodInd[i]; while ( ListArcPtr != NULL ) { if ( ListArcPtr->APtr[l]->Status == SET_FREE ) fprintf(OUTFILE, "(%2d,%2d,%2d) -> ( %f - %f * %d) = %10.6f\n", l+1, ListArcPtr->APtr[l]->i+1, ListArcPtr->APtr[l]->j+1, ListArcPtr->APtr[l]->xij, Ml[l], ListArcPtr->APtr[l]->yij, ListArcPtr->APtr[l]->Gamma_wij); ListArcPtr = ListArcPtr->Next; } } } fprintf(OUTFILE, "\n||Gamma||^2 = %f\n", SqrSum); */ return(SqrSum); } void First_Ck (int SizeOfProblem, int *Base, int *Displacement, int *k, float *Ck) { #if _STEP_ == _HO_SE_ *k = 1; *Ck = 2; #elif _STEP_ == _HE_WO_CR_ *Base = 2 * SizeOfProblem; *Displacement = 2 * SizeOfProblem; *k = 0; *Ck = 2; #endif } void Next_Ck (int *Base, int *Displacement, int *k, float *Ck) { (*k)++; #if _STEP_ == _HO_SE_ if ( (*k % 4) == 0 ) *Ck /= 2; if ( (*k % 100) == 0 ) *Ck = 2; #elif _STEP_ == _HE_WO_CR_ if ( *k >= *Base ) { *Ck /= 2; *Displacement /= 2; if ( *Displacement < THRESHOLD ) *Displacement = THRESHOLD; *Base += *Displacement; } #endif } void UpdateMultipliers (int n_Nodes, int n_Levels, ListArcPtrType *AdjaNodInd, ListNodPtrType *RootNodLev, float tk) { int i, l; ListArcPtrType ListArcPtr; ListNodPtrType ListNodPtr; for ( l = 0; l < n_Levels; l++ ) { ListNodPtr = RootNodLev[l]; while ( ListNodPtr != NULL) { if ( ListNodPtr->NPtr->Status == SET_FREE ) { ListNodPtr->NPtr->vi += tk * ListNodPtr->NPtr->Gamma_vi; if ( ListNodPtr->NPtr->vi < 0 ) ListNodPtr->NPtr->vi = 0; } ListNodPtr = ListNodPtr->Next; } for ( i = 0; i < n_Nodes; i++ ) { ListArcPtr = AdjaNodInd[i]; while ( ListArcPtr != NULL ) { if ( ListArcPtr->APtr[l]->Status == SET_FREE ) { ListArcPtr->APtr[l]->wij += tk * ListArcPtr->APtr[l]->Gamma_wij; if ( ListArcPtr->APtr[l]->wij < 0 ) ListArcPtr->APtr[l]->wij = 0; } ListArcPtr = ListArcPtr->Next; } } } }