Greedy Algorithm 和 Hill Algorithm解决Capacitated Facility Location Problem

Capacitated Facility Location Problem：

Suppose there are n facilities and m customers. We wish to choose:

（1） which of the n facilities to open

（2）the assignment of customers to facilities

（3）The objective is to minimize the sum of the opening cost 

（4）The total demand assigned to a facility must not exceed its capacity.

问题解决方法描述：

Capacitated Facility Location Problem (CFLP) is generalization of the Simple Plant Location Problem. In contrast to that we now suppose each facility can product limited quantity of produce. It is importance and quite prevalent assumption. And although mathematical models of those problems not  very differ, but solving methods for CFLP are more difficult. The most efficient  methods for CFPL are Lagrangen relaxation methods and the matrix column generation method.

Now we discrabe the mathematical model as integer programming problem. Let a set I = {1,..., I} give potential facility locations by production some uniform product. The number ci 0is the opening cost of facility at location  iI ,  Vi0 is the utmost value of prodaction at this location.

A set J ={1,…,J} assign clients that require service. For each pair (ij)  gij0 is the production and transportation costs and  pij0 is a value of product from fasility i needed to client j.

 Let us define the following notations:

Then the Capacitated Facility Location Problem may be written ：

1、Greedy Algorithm：

一、算法思想

贪心法的基本思路：
——从问题的某一个初始解出发逐步逼近给定的目标，以尽可能快的地求得更好的解。当达到某算法中的某一步不能再继续前进时，算法停止。
该算法存在问题：
1. 不能保证求得的最后解是最佳的；
2. 不能用来求最大或最小解问题；
3. 只能求满足某些约束条件的可行解的范围。
4. greedy optimal control应该就是用这种算法求得的最优化控制。

在该情景中的应用：

1、遍历所有customer，对于每一个customer，将该customer安排到所有facility的cost进行从小到大排序

2、对排序好的cost进行遍历，假如该facility剩余的容量满足customer的要求，则将该customer安排到该facility

实现code：

Greedy.h:

#ifndef GREEDY_H
#define GREEDY_H

#include <vector>
using namespace std;  

class Greedy {
public:
   	Greedy(int facilityNum, int customerNum, vector<double> capacity, vector<double> cost, vector<double> demand, vector<vector<double>> assignment);
   	~Greedy();
   	double getResult();
   	int* get_assignment_list();
	void showResult();
private:
	double getCost(int facilityNum, int customerNum);
	int facilityNum;
	int customerNum;
	vector<double> capacity;
	vector<double> cost;
	vector<double> demand;
	vector<vector<double>> assignment;
	int *open_status;	
	int *assignment_list;
	double *remain_capacity;

};
#endif

Greedy.cpp：

#include "Greedy.h"
#include <iostream>
#include <cstring>
#include <limits.h>
using namespace std;

Greedy::Greedy(int facilityNum, int customerNum, vector<double> capacity, vector<double> cost,
						 vector<double> demand, vector<vector<double>> assignment) {
	this -> facilityNum = facilityNum;
	this -> customerNum = customerNum;
	this -> capacity = capacity;
	this -> cost = cost;
	this -> demand = demand;
	this -> assignment = assignment;
	this -> open_status = new int[facilityNum];
	this -> assignment_list = new int[customerNum];
	this -> remain_capacity = new double[facilityNum];
	for(int i = 0; i < facilityNum; i++) {
		open_status[i] = 0;
		remain_capacity[i] = capacity[i];
	}
	for(int i = 0; i < customerNum; i++)
		assignment_list[i] = 0;
}

Greedy::~Greedy() {
	delete open_status;
	delete assignment_list;
	delete remain_capacity;
}

double Greedy::getResult() {

	//int open_status[facilityNum] = {0};
	//int assignment_list[customerNum] = {0};
	//double remain_capacity[facilityNum];
	//for(int i = 0; i < capacity.size(); i++)
	//	remain_capacity[i] = capacity[i];
	
	// 按照贪心算法预估的cost值，选取最小的cost去安排设备
	for(int count = 0; count < customerNum; count++) {
		// 获取最高性价比，评估cost最小的分配
		double leastCost = INT_MAX;
		int fac = -1;
		int cus = -1;
		for(int j = 0; j < facilityNum; j++) {
			for(int i = 0; i < customerNum; i++) {
				// customer 已经分配
				if(assignment_list[i] != 0) continue;
				// 设备容量不足
				if(demand[i] > remain_capacity[j]) continue;
				// 有容量，获取贪心算法评估代价Cost
				double cost_j_i = getCost(j, i);
				//cout << j+1 << " " << i+1 << " " << cost_j_i << endl; 
				if(cost_j_i < leastCost) {
					leastCost = cost_j_i;
					fac = j+1;
					cus = i+1;
				}
			}
		}
		// 没有可进行的分配
		if(fac == -1 || cus == -1 || leastCost == INT_MAX) {
			cout << "[ERROR]" << endl;
		}
		else {
			// 进行分配
			open_status[fac-1] = 1;
			remain_capacity[fac-1] -= demand[cus-1];
			assignment_list[cus-1] = fac;
		}
	}
	// 计算总共费用
	double totalCost = 0;
	// 设备启动费用
	for(int j = 0; j < facilityNum; j++) {
		//printf("%d ", open_status[j]);
		if(open_status[j] == 1) {
			totalCost += cost[j];
		}
		//cout << remain_capacity[j] << " ";
	}
	//cout << endl;

	//printf("\n");
	// customer 费用
	for(int i = 0; i < customerNum; i++) {
		//printf("%d ", assignment_list[i]);
		totalCost += assignment[assignment_list[i]-1][i];
	}
	//printf("\n");
	
	return totalCost;
}

// 返回贪心算法中评估Cost的值
double Greedy::getCost(int facilityNum, int customerNum) {
	// 开启设备分摊平均占用费用
	double facilityCost = cost[facilityNum] * demand[customerNum] / capacity[facilityNum];
	// 消费者使用费用
	double assignmentCost = assignment[facilityNum][customerNum];
	// 返回总评估费用
	return facilityCost + assignmentCost;
}


int* Greedy::get_assignment_list() {
	return assignment_list;
}

void Greedy::showResult() {
	for(int j = 0; j < facilityNum; j++) {
		if(remain_capacity[j] != capacity[j]) {
			printf("1 ");
		}
		else printf("0 ");
	}
	printf("\n");
	/*for(int j = 0; j < facilityNum; j++) {
		printf("%0.lf ", remain_capacity[j]);
	}
	printf("\n");*/
	// customer 费用
	for(int i = 0; i < customerNum; i++) {
		printf("%d ", assignment_list[i]);
	}
	printf("\n");
}

main.cpp:

#include <fstream>  
#include <string>
#include <iostream>
#include <sstream>
#include <vector>
#include "Greedy.h"
#include "Greedy.cpp"
#include <time.h>
using namespace std;  
  
int main() {
	cout << "BEGIN" << endl;
	// 计数读取71个数据集
	int count = 1;
	// 文件路径
	string baseFileName = "../Instances/";
	

	// 循环读取71个数据集
	while(count <= 71) {
		// 获取文件名
		stringstream sstream;
		sstream << count;
		string pNum;
		sstream >> pNum;
		string fileName = "p" + pNum;

		// 流读取文件
	    ifstream in(baseFileName+fileName);
	    ofstream out("../result/result.txt");

	  	// 有该文件
	    if(in) {
	    	vector<double> sample;
	    	double d;
	    	while(in >> d) {
	    		sample.push_back(d);
	    	}
	    	
	    	
	    	// fac 和 cus 数量数据
	    	int facilityNum = sample[0];
	    	int customerNum = sample[1];

	    	// fac 的 capacity 和 cost 数据
	    	vector<double> capacity;
	    	vector<double> cost;
	    	int begin = 2;
	    	for(int i = 0; i < facilityNum; i++) {
	    		capacity.push_back(sample[2*i+begin]);
	    		cost.push_back(sample[2*i+begin+1]);
	    	}

	    	// cus 的 demand 数据
	    	vector<double> demand;
	    	begin = 2+facilityNum*2;
	    	for(int i = 0; i < customerNum; i++) 
	    		demand.push_back(sample[i+begin]);

	    	// cus 的 assignment 数据
	    	vector<vector<double>> assignment;
	    	begin = 2+facilityNum*2+customerNum;
	    	for(int j = 0; j < facilityNum; j++) {
	    		vector<double> temp;
	    		for(int i = 0; i < customerNum; i++)
	    			temp.push_back(sample[j*customerNum+i+begin]);
	    		assignment.push_back(temp);
	    	}
	    	// 方法类
			printf("------Solution1------%s------fac:%d------cus:%d---\n", fileName.c_str(), facilityNum, customerNum);
			Solution1 s1 = Solution1(facilityNum, customerNum, capacity, cost, demand, assignment);	
			clock_t start,stop;
    		start = clock();
	    	double result = s1.getResult();
	    	s1.showResult();
	    	stop = clock();
			double dur = (double)stop-start;
			printf("[Time]---%.0lfms\n", dur);
	    	printf("[Result]---%lf\n", result);
	    }
	    // 没有该文件
	    else cout <<"No Test Data Files" << endl;
	    count++;
	}

	cout << "END" << endl;
	
    return 0;
}

Result table(给出部分结果)：

result table

instances	results	times(ms)
P1	9343	15
P2	7933	4
P3	10108	8
P15	13187	8
P16	17187	23
P33	15852	29

P1

P2

P3

P15

P16

P33

etc..

2、Hill Algorithm：

 其实就是the matrix column generation method.

算法解释：
从当前的节点开始，和周围的邻居节点的值进行比较。 如果当前节点是最大的，那么返回当前节点，作为最大值 (既山峰最高点)；反之就用最高的邻居节点来，替换当前节点，从而实现向山峰的高处攀爬的目的。如此循环直到达到最高点。

算法优缺点
优点
避免遍历，通过启发选择部分节点，从而达到提高效率的目的。

缺点

因为不是全面搜索，所以结果可能不是最佳。

算法伪代码：

function HILL-CLIMBING(problem) returns a state that is a local maximum
inputs: problem, a problem
  local variables: current, a node
                         neighbor, a node
  current <- MAKE-NODE(INITIAL-STATE[problem])
  loop do
neighbor <- a highest-valued successor of current
if VALUE[neighbor]<= VALUE[current] then return STATE[current]
current <- neighbor

 

实现code：

main.cpp
 

#include "individual.hpp"
#include <iostream>
#include <ctime>
#include <string>
#include <vector>
#include <fstream>
#include <cstdlib>
#include <cstdio>

using namespace std;

int main() {
  // Log file stream initializing //
  ofstream csvOutput;
  ofstream detailOutput;
  csvOutput.open("result.csv");
  detailOutput.open("detail.txt");
  csvOutput << ",Result,Time(s)\n";

  // Timer //
  time_t start;
  // Status monitor //
  int count = 0, scoreRecord;
  int* facilities = NULL;
  // Case info //
  int consumerNum, facilityNum;
  // Directory path //
  string rootDir = "../Instances/p";
  // Load each file and execute algorithm //
  for (int number = 1; number < 72; ++number) {
  
    try {
      Individual::init(rootDir + to_string(number));
    }
    catch(string str) {
      cout << str;
      return 0;
    }
    catch(exception e) {
      cout << "Error while initializing: " << e.what() << endl;
      return 0;
    }
    // Initializing for case //
    count = 0;
    Individual *bestScore = new Individual;
    consumerNum = Individual::getCustomerNum();
    facilityNum = Individual::getFacilityNum();

    if (facilities != nullptr) delete[] facilities;
    facilities = new int[facilityNum];
    for (int i = 0; i < facilityNum; ++i) facilities[i] = 0;
    // start timer //
    start = time(nullptr);

    while ((scoreRecord = bestScore->getCost()) == INT32_MAX) {
      delete bestScore;
      bestScore = new Individual();
      cout << "*";
    }


    cout << "*";
    do {
      auto temp3 = bestScore->getNeighborsWithTryMin();
      
      for (int i = 0; i < consumerNum; ++i) {
        if (bestScore->getCost() > temp3[i]->getCost()) {
          bestScore = temp3[i];
          bestScore->getCost();
        }
        else {
          delete temp3[i];
        }
      }
      delete []temp3;
      // Are the best still him? //
      if (bestScore->getCost() == scoreRecord) {
        time_t timeSpent = time(nullptr) - start;
        csvOutput << "p" << number << "," << scoreRecord
          << "," <<  timeSpent<< "\n";
        detailOutput << "Testcase: p" << number 
          << " Final Cost: " << scoreRecord 
          << " Time Spent: " << timeSpent << " second(s)" << endl;
        int *solution = bestScore->getSolution();
        for (int t = 0; t < consumerNum; ++t) {
          facilities[solution[t]] = 1;
        }

        for (int t = 0; t < facilityNum; ++t) {
          detailOutput << facilities[t] << " ";
        }
        detailOutput << endl;
        for (int t = 0; t < consumerNum; ++t) {
          detailOutput << solution[t] << " ";
        }
        detailOutput << endl;
        cout << "p" << number << " Done" << endl;
        break;
        
      } else {
        // Oh! A new king! //
        scoreRecord = bestScore->getCost();
        cout << "New Best At " << count << " :" << scoreRecord << endl;
        count = 0;
      }

    } while (true);
  }
  if (facilities != nullptr) delete[] facilities;
  detailOutput.close();
  csvOutput.close();
  return 0;
}

individual.h

#ifndef INDIVIDUAL_HPP
#define INDIVIDUAL_HPP
#include <string>

using std::string;


class Individual {
private:
  // Static variables from data file //
  static int    FACILITY_NUM;     // Number of facilities
  static int    CUSTOMER_NUM;     // Number of customers
  static int**  ASSIGN_COST;      // Matrix of customer's demand
  static int*   MIN_ASSIGN_COST;  // Minimal assign cost for each consumer
  static int*   MIN_ASSIGN;       // Minimal facility index
  static int*   CUSTOMER_DEMAND;  // Vector of facilities' assignment cost
  static int*   OPEN_COST;        // The cost of facility to open
  static int*   FACILITY_CAPACITY;// The capacity of each facility

  // Object's private data member //
  int*          gene;             // Gene, the strategy of assignment
  int           cost;

public:
  // Static methods //
  static void init(string filePath);              // Use data file to initialize static values
  static int  estimateCost(Individual &individual); // estimate grade for individual
  static void clearData();                        // Clear data and free memory space
  static bool isValid(int* gene);                 // Check whether a gene is valid

  static int getFacilityNum();
  static int getCustomerNum();

  // Constructors //
  Individual();
  ~Individual();
  Individual(int* gene);
  Individual(const Individual& individual);

  // Public methods //
  int*    getSolution();  // Get the solution
  int*    getGeneCopy();  // Get a copy of gene
  int     getCost();      // Get the cost
  Individual** getNeighborsWithTryMin();
  Individual& operator=(const Individual& individual);
};
#endif

individual.hpp

#include <fstream>
#include <cstdlib>
#include <algorithm>
#include "individual.hpp"

using std::ifstream;
using std::rand;
using std::swap;
#ifdef DEBUG_MODE
#include <iostream>
using std::cout;
using std::endl;
#endif

// Initialize static member variables in Individual //
int    Individual::FACILITY_NUM      = 0;     
int    Individual::CUSTOMER_NUM      = 0;     
int**  Individual::ASSIGN_COST       = nullptr;  
int*   Individual::CUSTOMER_DEMAND   = nullptr;      
int*   Individual::OPEN_COST         = nullptr; 
int*   Individual::FACILITY_CAPACITY = nullptr;
int*   Individual::MIN_ASSIGN_COST   = nullptr;
int*   Individual::MIN_ASSIGN        = nullptr;

// Clear all data and recycle memory space allocated //
void Individual::clearData() {
  if (ASSIGN_COST != nullptr){
    for (int i = 0; i < FACILITY_NUM; ++i) {
      delete[] ASSIGN_COST[i];
    }
    delete[] ASSIGN_COST;
    ASSIGN_COST = nullptr;
  }
  if (MIN_ASSIGN_COST != nullptr) {
    delete[] MIN_ASSIGN_COST;
    MIN_ASSIGN_COST = nullptr;
  }
  if (MIN_ASSIGN != nullptr) {
    delete[] MIN_ASSIGN;
    MIN_ASSIGN = nullptr;
  }
  if (CUSTOMER_DEMAND != nullptr) {
    delete[] CUSTOMER_DEMAND;
    CUSTOMER_DEMAND = nullptr;
  }
  if (OPEN_COST != nullptr) {
    delete[] OPEN_COST;
    OPEN_COST = nullptr;
  }
  if (FACILITY_CAPACITY != nullptr) {
    delete[] FACILITY_CAPACITY;
    FACILITY_CAPACITY = nullptr;
  }
}

// Initialize case info with the given file //
void Individual::init(string filePath) {
  // Create a input file stream //
  ifstream is;
  is.open(filePath); 
  double temp = 0;
  
  // If the file is not open properly //
  if (!is.is_open()) { 
    throw string("E: Unable To Open File");
    return;
  }

  // Clear existed data and free memory space
  clearData();
  
  // Number of facilities and customers //
  is >> FACILITY_NUM >> CUSTOMER_NUM;
  
  // Capacity and Open Cost //
  FACILITY_CAPACITY = new int[FACILITY_NUM];
  OPEN_COST = new int [FACILITY_NUM];
  for (int i = 0; i < FACILITY_NUM; ++i) {
    is >> FACILITY_CAPACITY[i] >> OPEN_COST[i];
  }

  // Demand of customers //
  CUSTOMER_DEMAND = new int[CUSTOMER_NUM];
  for (int i = 0; i < CUSTOMER_NUM; ++i) {
    is >> temp;
    CUSTOMER_DEMAND[i] = temp;
  }

  // Assign Cost //
  MIN_ASSIGN_COST = new int[CUSTOMER_NUM];
  MIN_ASSIGN = new int[CUSTOMER_NUM];
  for (int i = 0; i < CUSTOMER_NUM; ++i) MIN_ASSIGN_COST[i] = INT32_MAX;
  ASSIGN_COST = new int*[FACILITY_NUM];
  for (int i = 0; i < FACILITY_NUM; ++i) {
    ASSIGN_COST[i] = new int[CUSTOMER_NUM];
    for (int t = 0; t < CUSTOMER_NUM; ++t) {
      is >> temp;
      ASSIGN_COST[i][t] = temp;
      if (temp < MIN_ASSIGN_COST[t]){
        MIN_ASSIGN_COST[t] = temp;
        MIN_ASSIGN[t] = i;
      }
    }
  }

  // Debug info //
  #ifdef DEBUG_MODE
    cout << "Status after reading data:" << endl
      << "Facility Number: " << FACILITY_NUM << endl
      << "Customer Number: " << CUSTOMER_NUM << endl
      << "Assign Cost: " << endl;
    for (int i = 0; i < FACILITY_NUM; ++i) {
      for (int t = 0; t < CUSTOMER_NUM; ++t) {
        cout << ASSIGN_COST[i][t] << " ";
      }
      cout << endl;
    }
    cout << "Facility Capacities: " << endl;
    for (int i = 0; i < FACILITY_NUM; ++i) {
      cout << FACILITY_CAPACITY[i] << " ";
    }
    cout << endl
      << "Facility Open Cost: " << endl;
    for (int i = 0; i < FACILITY_NUM; ++i) {
      cout << OPEN_COST[i] << " ";
    } 
    cout << endl
      << "Customers' demand: " << endl;
    for (int i = 0; i < CUSTOMER_NUM; ++i) {
      cout << CUSTOMER_DEMAND[i] << " ";
    }
    cout << endl;
  #endif
}

// Estimate the cost of current solution //
int Individual::estimateCost(Individual& individual) {
  // If this solution is not a legal one, return the max cost //
  if (!isValid(individual.gene))
    return INT32_MAX;
  
  // Store the status of each facility //
  bool* facilities = new bool[FACILITY_NUM];
  int cost = 0;
  for (int i = 0; i < FACILITY_NUM; ++i) {
    facilities[i] = false;
  }

  // Count the cost //
  for (int i = 0; i < CUSTOMER_NUM; ++i) {
    if (facilities[individual.gene[i]] == false) {
      // Open the facility and count the cost //
      cost += OPEN_COST[individual.gene[i]];
      facilities[individual.gene[i]] = true;
    }
    // Count assignment cost //
    cost += ASSIGN_COST[individual.gene[i]][i];
  }

  // Free allocated memory space //
  delete []facilities;
  // Update and return the cost //
  return individual.cost = cost;
}

Individual::Individual() {
  // Brand new gene //
  this->gene = new int[CUSTOMER_NUM];
  
  // Sequence used to check the facility capacity //
  // int* facilities = new int[FACILITY_NUM];
  // for (int i = 0; i < FACILITY_NUM; ++i) {
  //   facilities[i] = 0;
  // }

  
  for (int i = 0; i < CUSTOMER_NUM; ++i) {
    // Random until valid //
    // do {
      // Assign a random facility to each customer //
      gene[i] = rand() % FACILITY_NUM;
    // } while (CUSTOMER_DEMAND[i] + 
    //   facilities[gene[i]] > FACILITY_CAPACITY[gene[i]]);
  }
  
  // Estimate the cost //
  cost = estimateCost(*this);
  // Free allocated memory space //
  // delete[] facilities;
}

Individual::~Individual() {
  // Free allocated memory space //
  if (gene != nullptr) delete[] gene;
}

int* Individual::getSolution() {
  return gene;
}

// Check whether this solution match requirement //
bool Individual::isValid(int* gene) {
  // Count demand //
  int* facilities = new int[FACILITY_NUM];
  bool flag = true;
  for (int i = 0; i < FACILITY_NUM; ++i) {
    facilities[i] = 0;
  }

  for (int i = 0; i < CUSTOMER_NUM; ++i) {
    facilities[gene[i]] += CUSTOMER_DEMAND[i];
    if (facilities[gene[i]] > FACILITY_CAPACITY[gene[i]]) {
      // Overflow! Not a valid solution //
      flag = false;
    }
  }

  // Free allocated memory //
  delete[] facilities;
  return flag;
}

Individual::Individual(const Individual& individual) {
  gene = new int[CUSTOMER_NUM];
  for (int i = 0; i < CUSTOMER_NUM; ++i) {
    // Copy gene //
    gene[i] = individual.gene[i];
  }
  // Copy cost //
  cost = individual.cost;
}

// Similar to copy constructor //
Individual& Individual::operator=(const Individual& individual) {
  if (gene == nullptr) gene = new int[CUSTOMER_NUM];
  for (int i = 0; i < CUSTOMER_NUM; ++i) {
    gene[i] = individual.gene[i];
  }
  cost = individual.cost;
  return *this;
}

// Get cost //
int Individual::getCost() {
  return cost;
}

// Get case info //
int Individual::getFacilityNum() {
  return FACILITY_NUM;
}
int Individual::getCustomerNum() {
  return CUSTOMER_NUM;
}

Individual** Individual::getNeighborsWithTryMin() {
  Individual** result = new Individual*[CUSTOMER_NUM];
  for (int i = 0; i < CUSTOMER_NUM; ++i) {
    auto temp = getGeneCopy();
    temp[i] = MIN_ASSIGN[i];
    result[i] = new Individual(temp);
  }
  return result;
}

// Get a copy of gene //
int* Individual::getGeneCopy() {
  int *copy = new int[CUSTOMER_NUM];
  for (int i = 0; i < CUSTOMER_NUM; ++i)
    copy[i] = gene[i];
  return copy;
}

Individual::Individual(int* gene) {
  // Copy gene //
  this->gene = gene;
  // Get Cost //
  cost = estimateCost(*this);
}

Result table(给出部分结果)：

result table

instances	results	times(ms)
P1	9546	24
P2	8890	36
p16	17922	56
p24	18322	11
P30	15001	21
P43	11913	44

P1

P2

P16

P24

p30

p43

etc..