allocation modelling problem

I am getting decentralized production from this code but I should get centralized production from one factory for the following parameter and cost function. I build the same model in excel and solved it and got centralized production. I am not sure what is wrong in my model. 

#from gurobipy import Model, GRB, quicksum
from gurobipy import Model, GRB, quicksum
import pandas as pd

# Define factories and markets
factories = [
   "H", "F1", "F2",
]

markets = [
   "H", "F1", "F2", "F3", 
]

# Parameters
G = {i: 10000 for i in factories}
a = {i: 500 for i in factories}
d = {
   "H": 1000,
   "F1": 1000,
   "F2": 1000,
   "F3": 1000,
}
k = {
   "H": 0.1,
   "F1":0.1,
   "F2": 0.1,
   "F3": 0.1,
}

# Model
model = Model("Profit_Maximization")
model.setParam("NonConvex", 2)

# Decision variables
z = model.addVars(factories, markets, lb=0, ub= 10000, name="z")

# Aggregate production per factory
s = model.addVars(factories, lb=0, name="s")       
s_pow = model.addVars(factories, lb=0, name="s_pow")

# Link factory total to shipments
for i in factories:
   model.addConstr(s[i] == quicksum(z[i, j] for j in markets), name=f"sum_{i}")

   alpha = -0.60
for i in factories:
   model.addGenConstrPow(s[i], s_pow[i], alpha, name=f"pow_{i}")

#  Revenue
Q = {j: d[j] - k[j] * quicksum(z[i, j] for i in factories) for j in markets}
Rtot = quicksum(Q[j] * quicksum(z[i, j] for i in factories) for j in markets)

# Cost 
# Factory-level economies of scale cost
C_power = quicksum(G[i] * s_pow[i] for i in factories)

# Per-lane costs (taxes, transport, etc.)
C_unit = quicksum(
   z[i, j] * a[i] 
   for i in factories for j in markets
)

Ctot = C_power + C_unit    

# Objective function
model.setObjective(Rtot - Ctot, GRB.MAXIMIZE)

# Solve
model.optimize()

# Calculate sales price per market
sales_prices = {}
total_quantities = {}

for j in markets:
   total_quantity = sum(z[i, j].X for i in factories)
   price = d[j] - k[j] * total_quantity
   sales_prices[j] = round(price, 2)
   total_quantities[j] = round(total_quantity, 2)

# Results
if model.status == GRB.OPTIMAL:
   # Total Revenue and Cost
   Rval = Rtot.getValue()
   Cval = Ctot.getValue()
   profit_margin = (Rval - Cval) / Rval if Rval != 0 else float('inf')

   print("\n✅ Maximum Total Profit (Ptot):", round(model.ObjVal, 2))
   print(f"📈 Total Revenue: {round(Rval, 2)}")
   print(f"💸 Total Cost: {round(Cval, 2)}")
   print(f"✅ Revenue based Profit Margin: {round(profit_margin * 100, 2)}%")
   print("\n🛒 Sales Price per Market and Supplied Quantity:")
   price_df = pd.DataFrame({
   "Sales Price": sales_prices,
   "Total Quantity Supplied": total_quantities
   })
   print(price_df)
   
   # Unit production cost per factory
   factory_costs = {}
   factory_quantities = {}
   unit_cost_per_factory = {}
   for i in factories:
       total_cost = 0
       # total production from this factory
       s_i = sum(z[i, j].X for j in markets)
       # factory-level costs
       econ_cost = G[i] * (s_i ** -0.60) if s_i > 0 else 0
       fixed_cost = a[i]
       # per-lane variable costs
       lane_cost = sum(
            z[i, j].X
            for j in markets
        )
       # total cost for this factory
       total_cost = econ_cost +fixed_cost+ lane_cost
       factory_costs[i] = round(total_cost, 2)
       # production quantity
       factory_quantities[i] = round(s_i, 2)
       # per-unit cost
       unit_cost_per_factory[i] = (
            round(total_cost / s_i, 2) if s_i > 0 else float("inf")
   )
   unit_cost_df = pd.DataFrame({
       "Total Cost": factory_costs,
       "Total Quantity": factory_quantities,
       "Unit Production Cost": unit_cost_per_factory
   })
   print("\n🏭 Unit Production Cost per Factory:")
   print(unit_cost_df)

   # Extract and show production matrix
   result = pd.DataFrame(index=markets, columns=factories)
   for i in factories:
       for j in markets:
           val = z[i, j].X
           result.loc[j, i] = round(val, 2) if val > 0.01 else 0.0

   print("\n📊 Production Matrix (units shipped from factories to markets):")
   print(result)
   

else:
   print("\n❌ No optimal solution found.")