From 8e987c3b77ecdce2ededc7f5ea522bbccd3d729c Mon Sep 17 00:00:00 2001 From: poveda_g Date: Thu, 11 Jan 2024 09:34:29 +0100 Subject: [PATCH 1/3] Adding modelisation of rcpsp using ortools library and cpsat solver. added associated test --- .../rcpsp/solver/cpsat_solver.py | 174 ++++++++++++++++++ tests/rcpsp/solver/test_rcpsp_cp.py | 27 +++ 2 files changed, 201 insertions(+) create mode 100644 discrete_optimization/rcpsp/solver/cpsat_solver.py diff --git a/discrete_optimization/rcpsp/solver/cpsat_solver.py b/discrete_optimization/rcpsp/solver/cpsat_solver.py new file mode 100644 index 000000000..a6c8d49ec --- /dev/null +++ b/discrete_optimization/rcpsp/solver/cpsat_solver.py @@ -0,0 +1,174 @@ +# Copyright (c) 2024 AIRBUS and its affiliates. +# This source code is licensed under the MIT license found in the +# LICENSE file in the root directory of this source tree. +# Native ortools-cpsat implementation for multimode rcpsp with resource calendar. +# Note : Model could most likely be improved with https://github.com/google/or-tools/blob/stable/examples/python/rcpsp_sat.py +import logging +from typing import Optional, Dict, Any, List +from discrete_optimization.generic_tools.result_storage.result_storage import ResultStorage, \ + from_solutions_to_result_storage +from discrete_optimization.rcpsp.rcpsp_model import RCPSPModel, RCPSPSolution +from discrete_optimization.generic_tools.do_solver import SolverDO +from discrete_optimization.generic_tools.cp_tools import ParametersCP, StatusSolver +from discrete_optimization.generic_tools.do_problem import (ParamsObjectiveFunction, + build_aggreg_function_and_params_objective) +from ortools.sat.python.cp_model import CpModel, CpSolver, VarArraySolutionPrinter, IntervalVar, \ + VarArrayAndObjectiveSolutionPrinter +from ortools.sat.python.cp_model import UNKNOWN, OPTIMAL, INFEASIBLE, FEASIBLE + +from discrete_optimization.rcpsp.rcpsp_utils import create_fake_tasks + +logger = logging.getLogger(__name__) + + +class CPSatRCPSPSolver(SolverDO): + def __init__(self, problem: RCPSPModel, + params_objective_function: Optional[ParamsObjectiveFunction] = None): + self.problem = problem + ( + self.aggreg_sol, + self.aggreg_from_dict_values, + self.params_objective_function, + ) = build_aggreg_function_and_params_objective( + self.problem, params_objective_function=params_objective_function + ) + self.cp_model: Optional[CpModel] = None + self.cp_solver: Optional[CpSolver] = None + self.variables: Optional[Dict[str, Any]] = None + self.status_solver: Optional[StatusSolver] = None + + def init_model(self): + model = CpModel() + starts_var = {} + ends_var = {} + is_present_var = {} + interval_var = {} + for task in self.problem.tasks_list: + starts_var[task] = model.NewIntVar(lb=0, ub=self.problem.horizon, + name=f"start_{task}") + ends_var[task] = model.NewIntVar(lb=0, ub=self.problem.horizon, + name=f"end_{task}") + interval_per_tasks = {} + for task in self.problem.mode_details: + interval_per_tasks[task] = set() + for mode in self.problem.mode_details[task]: + is_present_var[(task, mode)] = model.NewBoolVar(f"is_present_{task,mode}") + interval_var[(task, mode)] = \ + model.NewOptionalIntervalVar(start=starts_var[task], + size=self.problem.mode_details[task][mode]["duration"], + end=ends_var[task], + is_present=is_present_var[(task, mode)], + name=f"interval_{task,mode}") + interval_per_tasks[task].add((task, mode)) + # Precedence constraints + for task in self.problem.successors: + for successor_task in self.problem.successors[task]: + model.Add(starts_var[successor_task] >= ends_var[task]) + # 1 mode selected + for task in interval_per_tasks: + model.AddExactlyOne([is_present_var[k] for k in interval_per_tasks[task]]) + + resources = self.problem.resources_list + if self.problem.is_calendar: + fake_task: List[Dict[str, int]] = create_fake_tasks(rcpsp_problem=self.problem) + else: + fake_task = [] + for resource in resources: + if resource in self.problem.non_renewable_resources: + task_modes_consuming = [((task, mode), self.problem.mode_details[task][mode].get(resource, 0)) + for task in self.problem.tasks_list + for mode in self.problem.mode_details[task] + if self.problem.mode_details[task][mode].get(resource, 0) > 0] + model.Add(sum([is_present_var[x[0]]*x[1] + for x in task_modes_consuming]) <= self.problem.get_max_resource_capacity(resource)) + + else: + task_modes_consuming = [((task, mode), self.problem.mode_details[task][mode].get(resource, 0)) + for task in self.problem.tasks_list + for mode in self.problem.mode_details[task] + if self.problem.mode_details[task][mode].get(resource, 0) > 0] + fake_task_res = [(model.NewFixedSizeIntervalVar(start=f["start"], + size=f["duration"], + name=f"res_"), + f.get(resource, 0)) + for f in fake_task if f.get(resource, 0) > 0] + capacity = self.problem.get_max_resource_capacity(resource) + if capacity>1: + model.AddCumulative([interval_var[x[0]] + for x in task_modes_consuming]+[x[0] for x in fake_task_res], + demands=[x[1] for x in task_modes_consuming]+[x[1] for x in fake_task_res], + capacity=self.problem.get_max_resource_capacity(resource)) + if capacity == 1: + model.AddNoOverlap([interval_var[x[0]] + for x in task_modes_consuming] + [x[0] for x in fake_task_res]) + + model.Minimize(starts_var[self.problem.sink_task]) + self.cp_model = model + self.variables = {"start": starts_var, "end": ends_var, + "is_present": is_present_var} + + def solve(self, parameters_cp: Optional[ParametersCP] = None, **kwargs: Any) -> ResultStorage: + if self.cp_model is None: + self.init_model() + if parameters_cp is None: + parameters_cp = ParametersCP.default() + solver = CpSolver() + solver.parameters.max_time_in_seconds = parameters_cp.time_limit + solver.parameters.num_workers = parameters_cp.nb_process + callback = VarArrayAndObjectiveSolutionPrinter( + variables= + list(self.variables["is_present"].values())+list(self.variables["start"].values()) + ) + status = solver.Solve(self.cp_model, callback) + if status == UNKNOWN: + logger.info("Unknown") + self.status_solver = StatusSolver.UNKNOWN + if status == INFEASIBLE: + logger.info("Inf") + self.status_solver = StatusSolver.UNSATISFIABLE + if status == OPTIMAL: + logger.info("Optimal") + self.status_solver = StatusSolver.OPTIMAL + if status == FEASIBLE: + logger.info("Feasible") + self.status_solver = StatusSolver.SATISFIED + logger.info( + f"Solver finished, status={solver.StatusName(status)}, objective = {solver.ObjectiveValue()}," + f"best obj bound = {solver.BestObjectiveBound()}" + ) + return self.retrieve_solution(solver=solver) + + def retrieve_solution(self, solver: CpSolver): + schedule = {} + modes_dict = {} + for task in self.variables["start"]: + schedule[task] = {"start_time": solver.Value(self.variables["start"][task]), + "end_time": solver.Value(self.variables["end"][task])} + for task,mode in self.variables["is_present"]: + if solver.Value(self.variables["is_present"][task,mode]): + modes_dict[task] = mode + sol = RCPSPSolution(problem=self.problem, + rcpsp_schedule=schedule, + rcpsp_modes=[modes_dict[t] for t in self.problem.tasks_list_non_dummy]) + return from_solutions_to_result_storage([sol], + problem=self.problem, + params_objective_function=self.params_objective_function) + + +def cpstatus_to_dostatus(status_from_cpsat) -> StatusSolver: + """ + + :param status_from_cpsat: either [UNKNOWN,INFEASIBLE,OPTIMAL,FEASIBLE] from ortools.cp api. + :return: Status + """ + if status_from_cpsat == UNKNOWN: + return StatusSolver.UNKNOWN + if status_from_cpsat == INFEASIBLE: + return StatusSolver.UNSATISFIABLE + if status_from_cpsat == OPTIMAL: + return StatusSolver.OPTIMAL + if status_from_cpsat == FEASIBLE: + return StatusSolver.SATISFIED + + + diff --git a/tests/rcpsp/solver/test_rcpsp_cp.py b/tests/rcpsp/solver/test_rcpsp_cp.py index db0d1a26c..0732b6cc8 100644 --- a/tests/rcpsp/solver/test_rcpsp_cp.py +++ b/tests/rcpsp/solver/test_rcpsp_cp.py @@ -26,6 +26,7 @@ CP_MRCPSP_MZN_NOBOOL, CP_RCPSP_MZN, ) +from discrete_optimization.rcpsp.solver.cpsat_solver import CPSatRCPSPSolver @pytest.mark.parametrize( @@ -59,6 +60,32 @@ def test_cp_sm(optimisation_level): plot_task_gantt(rcpsp_problem, solution) +@pytest.mark.parametrize( + "model", + ["j301_1.sm", "j1010_1.mm"], +) +def test_ortools(model): + files_available = get_data_available() + file = [f for f in files_available if model in f][0] + rcpsp_problem = parse_file(file) + solver = CPSatRCPSPSolver(problem=rcpsp_problem) + parameters_cp = ParametersCP.default() + parameters_cp.time_limit = 100 + parameters_cp.nr_solutions = 1 + result_storage = solver.solve(parameters_cp=parameters_cp) + solution, fit = result_storage.get_best_solution_fit() + solution_rebuilt = RCPSPSolution( + problem=rcpsp_problem, + rcpsp_permutation=solution.rcpsp_permutation, + rcpsp_modes=solution.rcpsp_modes, + ) + fit_2 = rcpsp_problem.evaluate(solution_rebuilt) + assert fit == -fit_2["makespan"] + assert rcpsp_problem.satisfy(solution) + rcpsp_problem.plot_ressource_view(solution) + plot_task_gantt(rcpsp_problem, solution) + + def test_cp_sm_intermediate_solution(): files_available = get_data_available() file = [f for f in files_available if "j1201_1.sm" in f][0] From 2a244351a2e56a9bbfd7b327abb53b0418104587 Mon Sep 17 00:00:00 2001 From: poveda_g Date: Thu, 11 Jan 2024 09:39:34 +0100 Subject: [PATCH 2/3] rerun linter --- .../rcpsp/solver/cpsat_solver.py | 185 +++++++++++------- 1 file changed, 119 insertions(+), 66 deletions(-) diff --git a/discrete_optimization/rcpsp/solver/cpsat_solver.py b/discrete_optimization/rcpsp/solver/cpsat_solver.py index a6c8d49ec..f5c9b1f85 100644 --- a/discrete_optimization/rcpsp/solver/cpsat_solver.py +++ b/discrete_optimization/rcpsp/solver/cpsat_solver.py @@ -4,26 +4,42 @@ # Native ortools-cpsat implementation for multimode rcpsp with resource calendar. # Note : Model could most likely be improved with https://github.com/google/or-tools/blob/stable/examples/python/rcpsp_sat.py import logging -from typing import Optional, Dict, Any, List -from discrete_optimization.generic_tools.result_storage.result_storage import ResultStorage, \ - from_solutions_to_result_storage -from discrete_optimization.rcpsp.rcpsp_model import RCPSPModel, RCPSPSolution -from discrete_optimization.generic_tools.do_solver import SolverDO -from discrete_optimization.generic_tools.cp_tools import ParametersCP, StatusSolver -from discrete_optimization.generic_tools.do_problem import (ParamsObjectiveFunction, - build_aggreg_function_and_params_objective) -from ortools.sat.python.cp_model import CpModel, CpSolver, VarArraySolutionPrinter, IntervalVar, \ - VarArrayAndObjectiveSolutionPrinter -from ortools.sat.python.cp_model import UNKNOWN, OPTIMAL, INFEASIBLE, FEASIBLE +from typing import Any, Dict, List, Optional + +from ortools.sat.python.cp_model import ( + FEASIBLE, + INFEASIBLE, + OPTIMAL, + UNKNOWN, + CpModel, + CpSolver, + IntervalVar, + VarArrayAndObjectiveSolutionPrinter, + VarArraySolutionPrinter, +) +from discrete_optimization.generic_tools.cp_tools import ParametersCP, StatusSolver +from discrete_optimization.generic_tools.do_problem import ( + ParamsObjectiveFunction, + build_aggreg_function_and_params_objective, +) +from discrete_optimization.generic_tools.do_solver import SolverDO +from discrete_optimization.generic_tools.result_storage.result_storage import ( + ResultStorage, + from_solutions_to_result_storage, +) +from discrete_optimization.rcpsp.rcpsp_model import RCPSPModel, RCPSPSolution from discrete_optimization.rcpsp.rcpsp_utils import create_fake_tasks logger = logging.getLogger(__name__) class CPSatRCPSPSolver(SolverDO): - def __init__(self, problem: RCPSPModel, - params_objective_function: Optional[ParamsObjectiveFunction] = None): + def __init__( + self, + problem: RCPSPModel, + params_objective_function: Optional[ParamsObjectiveFunction] = None, + ): self.problem = problem ( self.aggreg_sol, @@ -44,21 +60,26 @@ def init_model(self): is_present_var = {} interval_var = {} for task in self.problem.tasks_list: - starts_var[task] = model.NewIntVar(lb=0, ub=self.problem.horizon, - name=f"start_{task}") - ends_var[task] = model.NewIntVar(lb=0, ub=self.problem.horizon, - name=f"end_{task}") + starts_var[task] = model.NewIntVar( + lb=0, ub=self.problem.horizon, name=f"start_{task}" + ) + ends_var[task] = model.NewIntVar( + lb=0, ub=self.problem.horizon, name=f"end_{task}" + ) interval_per_tasks = {} for task in self.problem.mode_details: interval_per_tasks[task] = set() for mode in self.problem.mode_details[task]: - is_present_var[(task, mode)] = model.NewBoolVar(f"is_present_{task,mode}") - interval_var[(task, mode)] = \ - model.NewOptionalIntervalVar(start=starts_var[task], - size=self.problem.mode_details[task][mode]["duration"], - end=ends_var[task], - is_present=is_present_var[(task, mode)], - name=f"interval_{task,mode}") + is_present_var[(task, mode)] = model.NewBoolVar( + f"is_present_{task,mode}" + ) + interval_var[(task, mode)] = model.NewOptionalIntervalVar( + start=starts_var[task], + size=self.problem.mode_details[task][mode]["duration"], + end=ends_var[task], + is_present=is_present_var[(task, mode)], + name=f"interval_{task,mode}", + ) interval_per_tasks[task].add((task, mode)) # Precedence constraints for task in self.problem.successors: @@ -70,44 +91,73 @@ def init_model(self): resources = self.problem.resources_list if self.problem.is_calendar: - fake_task: List[Dict[str, int]] = create_fake_tasks(rcpsp_problem=self.problem) + fake_task: List[Dict[str, int]] = create_fake_tasks( + rcpsp_problem=self.problem + ) else: fake_task = [] for resource in resources: if resource in self.problem.non_renewable_resources: - task_modes_consuming = [((task, mode), self.problem.mode_details[task][mode].get(resource, 0)) - for task in self.problem.tasks_list - for mode in self.problem.mode_details[task] - if self.problem.mode_details[task][mode].get(resource, 0) > 0] - model.Add(sum([is_present_var[x[0]]*x[1] - for x in task_modes_consuming]) <= self.problem.get_max_resource_capacity(resource)) + task_modes_consuming = [ + ( + (task, mode), + self.problem.mode_details[task][mode].get(resource, 0), + ) + for task in self.problem.tasks_list + for mode in self.problem.mode_details[task] + if self.problem.mode_details[task][mode].get(resource, 0) > 0 + ] + model.Add( + sum([is_present_var[x[0]] * x[1] for x in task_modes_consuming]) + <= self.problem.get_max_resource_capacity(resource) + ) else: - task_modes_consuming = [((task, mode), self.problem.mode_details[task][mode].get(resource, 0)) - for task in self.problem.tasks_list - for mode in self.problem.mode_details[task] - if self.problem.mode_details[task][mode].get(resource, 0) > 0] - fake_task_res = [(model.NewFixedSizeIntervalVar(start=f["start"], - size=f["duration"], - name=f"res_"), - f.get(resource, 0)) - for f in fake_task if f.get(resource, 0) > 0] + task_modes_consuming = [ + ( + (task, mode), + self.problem.mode_details[task][mode].get(resource, 0), + ) + for task in self.problem.tasks_list + for mode in self.problem.mode_details[task] + if self.problem.mode_details[task][mode].get(resource, 0) > 0 + ] + fake_task_res = [ + ( + model.NewFixedSizeIntervalVar( + start=f["start"], size=f["duration"], name=f"res_" + ), + f.get(resource, 0), + ) + for f in fake_task + if f.get(resource, 0) > 0 + ] capacity = self.problem.get_max_resource_capacity(resource) - if capacity>1: - model.AddCumulative([interval_var[x[0]] - for x in task_modes_consuming]+[x[0] for x in fake_task_res], - demands=[x[1] for x in task_modes_consuming]+[x[1] for x in fake_task_res], - capacity=self.problem.get_max_resource_capacity(resource)) + if capacity > 1: + model.AddCumulative( + [interval_var[x[0]] for x in task_modes_consuming] + + [x[0] for x in fake_task_res], + demands=[x[1] for x in task_modes_consuming] + + [x[1] for x in fake_task_res], + capacity=self.problem.get_max_resource_capacity(resource), + ) if capacity == 1: - model.AddNoOverlap([interval_var[x[0]] - for x in task_modes_consuming] + [x[0] for x in fake_task_res]) + model.AddNoOverlap( + [interval_var[x[0]] for x in task_modes_consuming] + + [x[0] for x in fake_task_res] + ) model.Minimize(starts_var[self.problem.sink_task]) self.cp_model = model - self.variables = {"start": starts_var, "end": ends_var, - "is_present": is_present_var} - - def solve(self, parameters_cp: Optional[ParametersCP] = None, **kwargs: Any) -> ResultStorage: + self.variables = { + "start": starts_var, + "end": ends_var, + "is_present": is_present_var, + } + + def solve( + self, parameters_cp: Optional[ParametersCP] = None, **kwargs: Any + ) -> ResultStorage: if self.cp_model is None: self.init_model() if parameters_cp is None: @@ -116,8 +166,8 @@ def solve(self, parameters_cp: Optional[ParametersCP] = None, **kwargs: Any) -> solver.parameters.max_time_in_seconds = parameters_cp.time_limit solver.parameters.num_workers = parameters_cp.nb_process callback = VarArrayAndObjectiveSolutionPrinter( - variables= - list(self.variables["is_present"].values())+list(self.variables["start"].values()) + variables=list(self.variables["is_present"].values()) + + list(self.variables["start"].values()) ) status = solver.Solve(self.cp_model, callback) if status == UNKNOWN: @@ -142,17 +192,23 @@ def retrieve_solution(self, solver: CpSolver): schedule = {} modes_dict = {} for task in self.variables["start"]: - schedule[task] = {"start_time": solver.Value(self.variables["start"][task]), - "end_time": solver.Value(self.variables["end"][task])} - for task,mode in self.variables["is_present"]: - if solver.Value(self.variables["is_present"][task,mode]): + schedule[task] = { + "start_time": solver.Value(self.variables["start"][task]), + "end_time": solver.Value(self.variables["end"][task]), + } + for task, mode in self.variables["is_present"]: + if solver.Value(self.variables["is_present"][task, mode]): modes_dict[task] = mode - sol = RCPSPSolution(problem=self.problem, - rcpsp_schedule=schedule, - rcpsp_modes=[modes_dict[t] for t in self.problem.tasks_list_non_dummy]) - return from_solutions_to_result_storage([sol], - problem=self.problem, - params_objective_function=self.params_objective_function) + sol = RCPSPSolution( + problem=self.problem, + rcpsp_schedule=schedule, + rcpsp_modes=[modes_dict[t] for t in self.problem.tasks_list_non_dummy], + ) + return from_solutions_to_result_storage( + [sol], + problem=self.problem, + params_objective_function=self.params_objective_function, + ) def cpstatus_to_dostatus(status_from_cpsat) -> StatusSolver: @@ -169,6 +225,3 @@ def cpstatus_to_dostatus(status_from_cpsat) -> StatusSolver: return StatusSolver.OPTIMAL if status_from_cpsat == FEASIBLE: return StatusSolver.SATISFIED - - - From 074a97dab68e3b572ce4dcd848d4435de6bd0972 Mon Sep 17 00:00:00 2001 From: poveda_g Date: Thu, 11 Jan 2024 10:14:50 +0100 Subject: [PATCH 3/3] transform cpsat solve status into discreteopt solve status --- discrete_optimization/rcpsp/solver/cpsat_solver.py | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/discrete_optimization/rcpsp/solver/cpsat_solver.py b/discrete_optimization/rcpsp/solver/cpsat_solver.py index f5c9b1f85..1c355c564 100644 --- a/discrete_optimization/rcpsp/solver/cpsat_solver.py +++ b/discrete_optimization/rcpsp/solver/cpsat_solver.py @@ -170,18 +170,7 @@ def solve( + list(self.variables["start"].values()) ) status = solver.Solve(self.cp_model, callback) - if status == UNKNOWN: - logger.info("Unknown") - self.status_solver = StatusSolver.UNKNOWN - if status == INFEASIBLE: - logger.info("Inf") - self.status_solver = StatusSolver.UNSATISFIABLE - if status == OPTIMAL: - logger.info("Optimal") - self.status_solver = StatusSolver.OPTIMAL - if status == FEASIBLE: - logger.info("Feasible") - self.status_solver = StatusSolver.SATISFIED + self.status_solver = cpstatus_to_dostatus(status_from_cpsat=status) logger.info( f"Solver finished, status={solver.StatusName(status)}, objective = {solver.ObjectiveValue()}," f"best obj bound = {solver.BestObjectiveBound()}"