基于CPLEX的化工厂反应器扩容决策建模技术咨询

Got it, let's break down how to integrate the reactor investment decision into your existing model. The core challenge is linking the "add/no-add + model choice" decision to dependent factors like initial stock, reactor capacity constraints, and net profit (accounting for reactor costs). Here's a step-by-step implementation:

1. Add Decision Variables for Reactor Selection #

First, define variables to represent your choice: no new reactor, or one of the three models. We'll use boolean variables to enforce exactly one option is selected.

// Define reactor options: 0 = no new reactor, 1/2/3 = three different models
{int} ReactorOptions = {0, 1, 2, 3};
dvar boolean ChooseReactor[m in ReactorOptions];

// Constraint: Must select exactly one option (no reactor or one model)
subject to {
    sum(m in ReactorOptions) ChooseReactor[m] == 1;
}

2. Calculate Total Reactors & Link to Initial Stock #

Your initial stock depends on how many reactors are running. First, compute the total number of reactors, then define initial stock as a function of your selection.

// Total reactors: original 3 + 1 if we choose any model (0 adds 0)
dexpr int TotalReactors = 3 + sum(m in ReactorOptions: m != 0) ChooseReactor[m];

// Define initial stock values per scenario:
// - InitialStock_3: stock when using original 3 reactors
// - InitialStock_1/2/3: stock when adding model 1/2/3
int InitialStock_3[Resources] = ...; // Your original initial stock values
int InitialStock_1[Resources] = ...; // Stock for model 1 addition
int InitialStock_2[Resources] = ...; // Stock for model 2 addition
int InitialStock_3_new[Resources] = ...; // Stock for model 3 addition

// Dynamic initial stock based on reactor choice
dexpr int InitialStock[r in Resources] = 
    InitialStock_3[r] * ChooseReactor[0] +
    InitialStock_1[r] * ChooseReactor[1] +
    InitialStock_2[r] * ChooseReactor[2] +
    InitialStock_3_new[r] * ChooseReactor[3];

3. Update Objective Function to Account for Reactor Costs #

Your goal is still maximizing net profit, so subtract the cost of the selected reactor from your original profit calculation.

// Cost of each reactor option: 0 = no cost for no reactor
float ReactorCost[m in ReactorOptions] = [0, cost_model1, cost_model2, cost_model3];

// Updated objective: original profit minus reactor cost
dexpr float ObjFunction = 
    sum(r in nrenuableR) (R[r][MaxTime] - InitialStock[r]) * (Profit[r] - nRcosts[r])
    - sum(m in ReactorOptions) ReactorCost[m] * ChooseReactor[m];
maximize ObjFunction;

4. Fix Reactor Scheduling Constraints #

You mentioned each reactor can only do one task at a time—your original model was missing this critical constraint. We'll add assignment variables to link tasks to reactors, and enforce exclusivity based on the total number of active reactors.

// Max possible reactors is 4 (original 3 + 1 new one)
range Reactors = 1..4;
// Assign[k][t][s] = 1 if task k starts at time t on reactor s
dvar boolean Assign[k in Tasks][t in Time][s in Reactors];
// ReactorActive[s] = 1 if reactor s is in use (depends on TotalReactors)
dvar boolean ReactorActive[s in Reactors];

// Enforce total active reactors matches our selection
subject to {
    sum(s in Reactors) ReactorActive[s] == TotalReactors;
}

// Tasks can only be assigned to active reactors
forall(k in Tasks, t in Time, s in Reactors) {
    Assign[k][t][s] <= ReactorActive[s];
}

// A reactor can't run more than one task at the same time
forall(s in Reactors, t in Time) {
    sum(k in Tasks) Assign[k][t][s] <= 1;
}

// If a task starts at time t (N[k][t] = 1), it must be assigned to a reactor
forall(k in Tasks, t in Time) {
    N[k][t] <= sum(s in Reactors) Assign[k][t][s];
}

// Prevent overlapping tasks on the same reactor (account for process time)
forall(k in Tasks, t in Time, s in Reactors) {
    forall(tp in t..min(t + procTime[k] - 1, MaxTime)) {
        Assign[k][t][s] <= 1 - sum(k2 in Tasks: k2 != k) Assign[k2][tp][s];
    }
}

Full Modified Code #

Here's how all the pieces fit together (integrated with your original code):

// Data declaration
int MaxTime = ...;
range Time = 0..MaxTime;
{int} Tasks = ...;
{string} nrenuableR=...;
{string} renuableR=...;
{string} renuableRusedbyT[Tasks]=...;
{string} Resources= nrenuableR union renuableR;
int procTime[Tasks]= ...;
int minbatchsize[renuableR][Tasks] =...;
int maxbatchsize [renuableR][Tasks] =...;
int MaxAmountStock_nR[nrenuableR]=...;
int maxRenuableR[renuableR][Time] =...;
int Profit[nrenuableR]=...;
float nRcosts[nrenuableR]=...;
int MaxTheta = ...;
range Theta=0..MaxTheta;
float Mu[Tasks][Resources][Theta] = ...;
float Nu[Tasks][Resources][Theta] = ...;

// Reactor investment data
{int} ReactorOptions = {0, 1, 2, 3};
float ReactorCost[m in ReactorOptions] = [0, cost_model1, cost_model2, cost_model3];
int InitialStock_3[Resources] = ...; // Original 3 reactors stock
int InitialStock_1[Resources] = ...; // Stock with model 1 added
int InitialStock_2[Resources] = ...; // Stock with model 2 added
int InitialStock_3_new[Resources] = ...; // Stock with model 3 added

// Decision variables
dvar boolean N[Tasks][Time];
dvar float+ Csi[Tasks][Time];
dvar int+ R[Resources][Time];
dvar boolean ChooseReactor[m in ReactorOptions];
range Reactors = 1..4;
dvar boolean Assign[k in Tasks][t in Time][s in Reactors];
dvar boolean ReactorActive[s in Reactors];

// Derived expressions
dexpr int TotalReactors = 3 + sum(m in ReactorOptions: m != 0) ChooseReactor[m];
dexpr int InitialStock[r in Resources] = 
    InitialStock_3[r] * ChooseReactor[0] +
    InitialStock_1[r] * ChooseReactor[1] +
    InitialStock_2[r] * ChooseReactor[2] +
    InitialStock_3_new[r] * ChooseReactor[3];

// Objective function
dexpr float ObjFunction = 
    sum(r in nrenuableR) (R[r][MaxTime] - InitialStock[r]) * (Profit[r] - nRcosts[r])
    - sum(m in ReactorOptions) ReactorCost[m] * ChooseReactor[m];
maximize ObjFunction;

// Constraint conditions
subject to {
    // Reactor selection constraint
    sum(m in ReactorOptions) ChooseReactor[m] == 1;
    
    // Reactor activation count constraint
    sum(s in Reactors) ReactorActive[s] == TotalReactors;
    
    // Resource capacity constraints
    forall (r in renuableR) forall(t in Time) R[r][t] <= maxRenuableR[r][t];
    forall (r in nrenuableR) forall (t in Time) R[r][t] <= MaxAmountStock_nR[r];
    
    // Batch size constraints
    forall (k in Tasks, r in renuableRusedbyT[k], t in Time) minbatchsize[r][k] * N[k][t] <= Csi[k][t];
    forall (k in Tasks, r in renuableRusedbyT[k], t in Time) maxbatchsize[r][k]*N[k][t] >= Csi[k][t];
    
    // Resource balance constraints
    forall(r in Resources) R[r][0] == InitialStock[r] + sum(k in Tasks) (Mu[k][r][0] * N[k][0] + Nu[k][r][0] * Csi[k][0]);
    forall(r in Resources,t in Time: t>0) R[r][t] == R[r][t-1] + sum(k in Tasks,theta in Theta: t - theta >=0) (Mu[k][r][theta] * N[k][t - theta] + Nu[k][r][theta] * Csi[k][t - theta]);
    
    // Reactor assignment constraints
    forall(k in Tasks, t in Time, s in Reactors) {
        Assign[k][t][s] <= ReactorActive[s];
    }
    forall(s in Reactors, t in Time) {
        sum(k in Tasks) Assign[k][t][s] <= 1;
    }
    forall(k in Tasks, t in Time) {
        N[k][t] <= sum(s in Reactors) Assign[k][t][s];
    }
    forall(k in Tasks, t in Time, s in Reactors) {
        forall(tp in t..min(t + procTime[k] - 1, MaxTime)) {
            Assign[k][t][s] <= 1 - sum(k2 in Tasks: k2 != k) Assign[k2][tp][s];
        }
    }
}

Key Notes #

• The ChooseReactor variables let CPLEX evaluate all four scenarios (no reactor + three models) and pick the one that maximizes net profit.
• Dynamic initial stock ensures you're using the correct starting inventory for each reactor count scenario.
• The reactor scheduling constraints enforce your "one task per reactor at a time" rule, even when the total number of reactors changes.
• Subtracting reactor costs from the objective ensures you're comparing true net profit across all options.

内容的提问来源于stack exchange，提问作者Beatriz Estevens