Solve a MILP model with Simulated Annealing in Python

Question

I applied a MILP mathematical programming model that considers integer and continuous variables and all functions (objective and constraints) are linear in pulp (applied with success).

nt = 5 # Periods
ni =  2 # Consumers
deltat = 0.25

piBuy = {(0, 0): 0.1034, (0, 1): 0.1034, (1, 0): 0.1034, (1, 1): 0.1034, (2, 0): 0.1034, (2, 1): 0.1034, (3, 0): 0.1034, (3, 1): 0.1034, (4, 0): 0.1034, (4, 1): 0.1034 }
piSell = {(0, 0): 0.045, (0, 1): 0.045, (1, 0): 0.045, (1, 1): 0.045, (2, 0): 0.045, (2, 1): 0.045, (3, 0): 0.045, (3, 1): 0.045, (4, 0): 0.045, (4, 1): 0.045 }
PgenPV = {(0, 0): 0, (0, 1): 0, (1, 0): 0, (1, 1): 0, (2, 0): 0, (2, 1): 0, (3, 0): 0, (3, 1): 0, (4, 0): 0, (4, 1): 0 }
Pload = {(0, 0): 0.22864236717888, (0, 1): 3.912695574263758, (1, 0): 0.341368689043474, (1, 1): 0.279951617541324, (2, 0): 0.248263477430366, (2, 1): 0.248014847666063, (3, 0): 0.277896110866711, (3, 1): 0.300043899274291, (4, 0): 0.33425916227626, (4, 1): 0.192176996274082 }
cFix = {0: 0.4729, 1: 0.6249}
maxPbuy = {0: 10.35, 1: 13.8}
maxPsell = {0: 10.35, 1: 13.8}
maxPch = {0: 2, 1: 0}
maxPdch = {0: 2, 1: 0}
upper = 1000000000

## Model
m = LpProblem("problem", LpMinimize)

# Create Decision variables
pBuy = {(t, i): LpVariable(f"pBuyGrid[{t},{i}]", lowBound=0, cat=LpContinuous) for t in range(nt) for i in range(ni)}
pSell = {(t, i): LpVariable(f"pSellGrid[{t},{i}]", lowBound=0, cat=LpContinuous) for t in range(nt) for i in range(ni)}
xBuy = {(t, i): LpVariable(f"binBuyGrid[{t},{i}]", cat=LpBinary) for t in range(nt) for i in range(ni)}
xSell = {(t, i): LpVariable(f"binSellGrid[{t},{i}]", cat=LpBinary) for t in range(nt) for i in range(ni)}
pDch = {(t, i): LpVariable(f"pdch[{t},{i}]", lowBound=0, cat=LpContinuous) for t in range(nt) for i in range(ni)}
pCh = {(t, i): LpVariable(f"pch[{t},{i}]", lowBound=0, cat=LpContinuous) for t in range(nt) for i in range(ni)}
xDch = {(t, i): LpVariable(f"binDch[{t},{i}]", cat=LpBinary) for t in range(nt) for i in range(ni)}
xCh = {(t, i): LpVariable(f"binCh[{t},{i}]", cat=LpBinary) for t in range(nt) for i in range(ni)}
pSellNo = {(t, i): LpVariable(f"pSellNo[{t},{i}]", lowBound=0, cat=LpContinuous) for t in range(nt) for i in range(ni)}
pExp = {(t, i): LpVariable(f"pExp[{t},{i}]", lowBound=0, cat=LpContinuous) for t in range(nt) for i in range(ni)}

# Objective function
m += lpSum((piBuy[t, i] * pBuy[t, i] - piSell[t, i] * pSell[t, i]) * deltat for t in range(nt) for i in range(ni)) + lpSum(cFix[i] for i in range(ni))

# Add constraints
for t in range(nt):
   for i in range(ni):
# constraint 1
        m += PgenPV[t,i] + pBuy[t, i] + pDch[t, i] == Pload[t,i] + pExp[t, i] + pCh[t, i]
# constraint 2
        m += pExp[t, i] == pSell[t, i] + pSellNo[t, i]
# limits
        m += pBuy[t, i] >= 0
        m += pBuy[t, i] <= maxPbuy[i] * xBuy[t, i]       
        m += pSell[t, i] >= 0
        m += pSell[t, i] <= maxPsell[i] * xSell[t, i]        
        m += xBuy[t, i] + xSell[t, i] <= 1 
        m += pSellNo[t, i] >= 0
        m += pSellNo[t, i] <= upper * xSell[t, i]
# batteries limits
        m += pCh[t, i] >= 0
        m += pCh[t, i] <= maxPch[i] * xCh[t, i] 
        m += pDch[t, i] >= 0
        m += pDch[t, i] <= maxPdch[i] * xDch[t, i] 
        m += xCh[t, i] + xDch[t, i] <= 1
# Dimension
        m += xBuy[t, i] >= 0
        m += xBuy[t, i] <= 1
        m += xSell[t, i] >= 0
        m += xSell[t, i] <= 1
        m += xCh[t, i] >= 0
        m += xCh[t, i] <= 1
        m += xDch[t, i] >= 0
        m += xDch[t, i] <= 1

# # Solve the problem
m.solve()
print("Costs:", value(m.objective))

The result returns the optimal value of: Costs: 1.1289073611495

Now I was trying to migrate the problem to be solved with a metaheuristic, in this case simulated annealing. I think I managed to model the main part of the simulated annealing approach in Python, but I'm having difficulty generating "perturbations in my sampling" and ensuring that the restrictions are met. I know I may have to add penalty functions to my objective function, but at this stage I wanted to try to understand what I'm doing wrong.

nt = 5 # Periods
ni =  2 # Consumers
deltat = 0.25

piBuy = {(0, 0): 0.1034, (0, 1): 0.1034, (1, 0): 0.1034, (1, 1): 0.1034, (2, 0): 0.1034, (2, 1): 0.1034, (3, 0): 0.1034, (3, 1): 0.1034, (4, 0): 0.1034, (4, 1): 0.1034 }
piSell = {(0, 0): 0.045, (0, 1): 0.045, (1, 0): 0.045, (1, 1): 0.045, (2, 0): 0.045, (2, 1): 0.045, (3, 0): 0.045, (3, 1): 0.045, (4, 0): 0.045, (4, 1): 0.045 }
PgenPV = {(0, 0): 0, (0, 1): 0, (1, 0): 0, (1, 1): 0, (2, 0): 0, (2, 1): 0, (3, 0): 0, (3, 1): 0, (4, 0): 0, (4, 1): 0 }
Pload = {(0, 0): 0.22864236717888, (0, 1): 3.912695574263758, (1, 0): 0.341368689043474, (1, 1): 0.279951617541324, (2, 0): 0.248263477430366, (2, 1): 0.248014847666063, (3, 0): 0.277896110866711, (3, 1): 0.300043899274291, (4, 0): 0.33425916227626, (4, 1): 0.192176996274082 }
cFix = {0: 0.4729, 1: 0.6249}
maxPbuy = {0: 10.35, 1: 13.8}
maxPsell = {0: 10.35, 1: 13.8}
maxPch = {0: 2, 1: 0}
maxPdch = {0: 2, 1: 0}
upper = 1000000000

# Initialize the state variables pBuy, pSell and others with random values between.
def initialize_states(nt, ni, maxPbuy, maxPsell, upper, maxPch, maxPdch):
    pBuy = {(t, i): max(random.uniform(0, maxPbuy[i]), 0) for t in range(nt) for i in range(ni)}
    pSell = {(t, i): max(random.uniform(0, maxPsell[i]), 0) for t in range(nt) for i in range(ni)}
    xBuy = {(t, i): random.choice([0, 1]) for t in range(nt) for i in range(ni)}
    xSell = {(t, i): random.choice([0, 1]) for t in range(nt) for i in range(ni)}
    pDch = {(t, i): max(random.uniform(0, maxPdch[i]), 0) for t in range(nt) for i in range(ni)}
    pCh = {(t, i): max(random.uniform(0, maxPch[i]), 0) for t in range(nt) for i in range(ni)}
    xCh = {(t, i): random.choice([0, 1]) for t in range(nt) for i in range(ni)}
    xDch = {(t, i): random.choice([0, 1]) for t in range(nt) for i in range(ni)}
    pSellNo = {(t, i): max(random.uniform(0, upper), 0) for t in range(nt) for i in range(ni)}
    pExp = {(t, i): pSell[t, i] + pSellNo[t, i] for t in range(nt) for i in range(ni)}
    
    return pBuy, pSell, xBuy, xSell, pDch, pCh, pExp, xCh, xDch, pSellNo

# Calculate the objective function value based on pBuy and pSell.
def objective_function(pBuy, pSell, piBuy, piSell, deltat, cFix):
    cost = sum((piBuy[t, i] * pBuy[(t, i)] - piSell[t, i] * pSell[(t, i)]) * deltat for t in range(nt) for i in range(ni))
    cost += sum(cFix[i] for i in range(ni))

    return cost

# Perturb the current solution by adding random values within the specified range.
def perturb(pBuy, pSell, perturbation_scale, lower_bound, maxPbuy, maxPsell, xBuy, xSell, pDch, pCh, pExp, xCh, xDch, pSellNo):

    new_pBuy = pBuy.copy()
    new_pSell = pSell.copy()
    new_xBuy = xBuy.copy()
    new_xSell = xSell.copy()
    new_pDch = pDch.copy()
    new_pCh = pCh.copy()
    new_pExp = pExp.copy()
    new_xCh = xCh.copy()
    new_xDch = xDch.copy()
    new_pSellNo = pSellNo.copy()

    for t in range(nt):
        for i in range(ni):

            # Configuration xBuy e xSell
            current_xBuy = xBuy[(t, i)]
            current_xSell = xSell[(t, i)]
            current_xCh = xCh[(t, i)]
            current_xDch = xDch[(t, i)]

            if current_xBuy == 1 and current_xSell == 1:
                if random.random() < 0.5:
                    new_xBuy[(t, i)] = 0
                else:
                    new_xSell[(t, i)] = 0

            if current_xCh == 1 and current_xDch == 1:
                if random.random() < 0.5:
                    new_xCh[(t, i)] = 0
                else:
                    new_xDch[(t, i)] = 0

            perturbation = random.uniform(0, perturbation_scale)

            new_pBuy[(t, i)] += perturbation
            new_pSell[(t, i)] += perturbation
            new_pDch[(t, i)] += perturbation
            new_pCh[(t, i)] += perturbation
            new_pExp[(t, i)] += perturbation
            new_pSellNo[(t, i)] += perturbation

            upper_bound_pBuy = maxPbuy[i] * new_xBuy[(t, i)]
            upper_bound_pSell = maxPsell[i] * new_xSell[(t, i)]
            upper_bound_pCh = maxPch[i] * new_xCh[t, i]
            upper_bound_pDch = maxPdch[i] * new_xDch[t, i] 
            upper_bound_pSellNo = upper * new_xSell[t, i]

            new_pBuy[(t, i)] = max(new_pBuy[(t, i)], 0)
            new_pSell[(t, i)] = max(new_pSell[(t, i)], 0)
            new_pDch[(t, i)] = max(new_pDch[(t, i)], 0)
            new_pCh[(t, i)] = max(new_pCh[(t, i)], 0)
            new_pExp[(t, i)] = max(new_pExp[(t, i)], 0)
            new_pSellNo[(t, i)] = max(new_pSellNo[(t, i)], 0)

            new_pBuy[(t, i)] = min(new_pBuy[(t, i)], upper_bound_pBuy)
            new_pSell[(t, i)] = min(new_pSell[(t, i)], upper_bound_pSell)
            new_pCh[(t, i)] = min(new_pCh[(t, i)], upper_bound_pCh)
            new_pDch[(t, i)] = min(new_pDch[(t, i)], upper_bound_pDch)
            new_pSellNo[(t, i)] = min(new_pSellNo[(t, i)], upper_bound_pSellNo)
            new_pExp[(t, i)] = new_pSell[(t, i)] + new_pSellNo[(t, i)]
    
            # First constraint
            left_side_sum = PgenPV[t, i] + pBuy[t, i] + pDch[t, i]
            right_side_sum = Pload[t, i] + pExp[t, i] + pCh[t, i]

            if not math.isclose(left_side_sum, right_side_sum, rel_tol=1e-9):
                total_difference = right_side_sum - left_side_sum

                # distribute the difference
                num_variables = 4 
                difference_per_variable = total_difference / num_variables

                # Adjust
                new_pBuy[t, i] += difference_per_variable
                new_pDch[t, i] += difference_per_variable
                new_pExp[t, i] += difference_per_variable
                new_pCh[t, i] += difference_per_variable

    return new_pBuy, new_pSell, new_xBuy, new_xSell, new_pCh, new_pDch, new_pSellNo, new_pExp

# Determine whether to accept a new solution based on the temperature and the cost difference.
def accept(old_cost, new_cost, temperature, scale=1.0):
    if new_cost < old_cost:
        return True
    else:
        probability = math.exp((old_cost - new_cost) / (temperature * scale))
        return random.random() < probability

# Simulated Annealing algorithm to optimize the solution.
def simulated_annealing(initial_temperature, cooling_rate, max_iterations, perturbation_scale, nt, ni, piBuy, piSell, deltat, cFix):
    pBuy, pSell, xBuy, xSell, pDch, pCh, pExp, xCh, xDch, pSellNo = initialize_states(nt, ni, maxPbuy, maxPsell, upper, maxPch, maxPdch)
    current_cost = objective_function(pBuy, pSell, piBuy, piSell, deltat, cFix)
    best_pBuy = pBuy
    best_pSell = pSell
    best_cost = current_cost
    temperature = initial_temperature

    costs = []  # List to store the costs at each iteration
    lower_bound = 0

    for i in range(max_iterations):
        new_pBuy, new_pSell, new_xBuy, new_xSell, new_pCh, new_pDch, new_pSellNo, new_pExp = perturb(pBuy, pSell, perturbation_scale, lower_bound, maxPbuy, maxPsell, xBuy, xSell, pDch, pCh, pExp, xCh, xDch, pSellNo)
        
        new_cost = objective_function(new_pBuy, new_pSell, piBuy, piSell, deltat, cFix)

        if accept(current_cost, new_cost, temperature):
            pBuy = new_pBuy
            pSell = new_pSell
            current_cost = new_cost

            if new_cost < best_cost:
                best_pBuy = new_pBuy
                best_pSell = new_pSell
                best_cost = new_cost

        temperature *= cooling_rate

        costs.append(best_cost)  # Append the cost to the list

        # Plot the cost over iterations
    plt.plot(range(max_iterations), costs)
    plt.xlabel("Iteration")
    plt.ylabel("Cost")
    plt.title("Cost vs. Iteration")
    plt.show()

    return best_pBuy, best_pSell, best_cost

initial_temperature = 100.0
cooling_rate = 0.98
max_iterations = 100
perturbation_scale = 0.05

best_pBuy, best_pSell, best_cost = simulated_annealing(initial_temperature, cooling_rate, max_iterations, perturbation_scale, nt, ni, piBuy, piSell, deltat, cFix)

Any idea what I'm doing wrong? the result is larger and appears to have no changes since the first iteration, which is strange.

I think the function of initializing the state of each variable is ok, I use random values and take into account whether it is variable or continuous, and the limits associated with each decision variable. My biggest difficulty is dealing with the disturbance function, where I am generating the samplings and including the restrictions:

Problem Constraints

I was hoping to reach a value close to what I get in pulp and so far I haven't achieved it. Any advice is welcome, thanks.