公司含硫原料生产

问题

   某公司将4种不同含硫量的液体原料(分别记为甲、乙、丙、丁)混合生产两种产品
   (分别记为A,B).按照生产工艺的要求,原料甲,乙、丁必须首先倒入混合池中混合,混
   合后的液体再分别与原料丙混合生产A,B.已知原料甲、乙,丙、丁的含硫量分别是
   3% , 1 % ,2% , 1 % ,进货价格(单位:千元/t)分别为6,16,10,15;产品A,B的含硫量分
   别不能超过2.5% ,1.5% ,售价分别为9,15.根据市场信息,原料甲、乙,丙的供应没有
   限制,原料丁的供应量最多为50t;产品A,B的市场需求量分别为100 t ,200 t.问应如何
   安排生产？

假设$y_1,z_1$分别为来自混合池和丙池的生成A的原料质量，$y_2,z_2$分别为来自混合池和丙池的生成B的原料质量，混合池中$x_1,x_2,x_4$分别为甲乙丁所占的比例。
那么可构建模型如下

$$\begin{gathered} Max \\ qquady=9\ast (y_1+z_1)+15\ast (y_2+z_2)-10\ast (z_1+z_2)-6\ast x_1\ast (y_1+y_2)-16\ast x_2\ast (y_1+y_2)-15\ast x_4\ast (y_1+y_2) \\ S.T.=\{\begin{array}{r}{x}_4\ast (y_1+y_2)<=50\\ y_1+z_1<=100\\ y_2+z_2<=200\\ (3\ast x_1+1\ast x_2+1\ast x_4)\ast y_1+2\ast z_1<=2.5\ast (y_1+z_1)\\ (3\ast x_1+1\ast x_2+1\ast x_4)\ast y_2+2\ast z_2<=1.5\ast (y_2+z_2)\\ x_1+x_2+x_4==1\end{array} \end{gathered}$$

import gurobipy as grb
from gurobipy import *
model = grb.Model("Make")
# 定义变量
y1 = model.addVar(name='y1')
y2 = model.addVar(name='y2')
z1 = model.addVar(name='z1')
z2 = model.addVar(name='z2')
x1 = model.addVar(name='x1')
x2 = model.addVar(name='x2')
x4 = model.addVar(name='x4')
# 添加约束
model.addConstr(x4*(y1+y2)<=50,"丁的供应量")
model.addConstr(y1+z1<=100,"A的需求量")
model.addConstr(y2+z2<=200,"B的需求量")
model.addConstr((3*x1+1*x2+1*x4)*y1+2*z1<=2.5*(y1+z1),"A含S占比")
model.addConstr((3*x1+1*x2+1*x4)*y2+2*z2<=1.5*(y2+z2),"B含S占比")
model.addConstr(x1+x2+x4==1,"比例约束")
model.Params.NonConvex = 2#将线性规划变为非线性规划 模型将按照求解MIP的逻辑来求解此类非凸问题。
model.setObjective( 9*(y1+z1)+15*(y2+z2)-10*(z1+z2)-6*x1*(y1+y2)-16*x2*(y1+y2)-15*x4*(y1+y2), grb.GRB.MAXIMIZE)
model.optimize()
if model.status == GRB.OPTIMAL:
    model.printAttr('X')
    model.printAttr('Slack')

Changed value of parameter NonConvex to 2
   Prev: -1  Min: -1  Max: 2  Default: -1
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 3 rows, 7 columns and 7 nonzeros
Model fingerprint: 0x34e4e21e
Model has 6 quadratic objective terms
Model has 3 quadratic constraints
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [1e+00, 3e+00]
  QLMatrix range   [5e-01, 3e+00]
  Objective range  [1e+00, 2e+01]
  QObjective range [1e+01, 3e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 2e+02]
  QRHS range       [5e+01, 5e+01]

Continuous model is non-convex -- solving as a MIP.

Presolve time: 0.00s
Presolved: 19 rows, 14 columns, 50 nonzeros
Presolved model has 6 bilinear constraint(s)
Variable types: 14 continuous, 0 integer (0 binary)

Root relaxation: objective 2.933333e+03, 7 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 2933.33333    0    6          - 2933.33333      -     -    0s
     0     0  550.00000    0    2          -  550.00000      -     -    0s
     0     0  550.00000    0    2          -  550.00000      -     -    0s
     0     0  525.00000    0    3          -  525.00000      -     -    0s
     0     0  450.00000    0    3          -  450.00000      -     -    0s
H    0     0                     450.0000000  450.00000  0.00%     -    0s

Cutting planes:
  RLT: 6

Explored 1 nodes (19 simplex iterations) in 0.09 seconds
Thread count was 8 (of 8 available processors)

Solution count 1: 450 

Optimal solution found (tolerance 1.00e-04)
Best objective 4.500000000000e+02, best bound 4.500000000000e+02, gap 0.0000%

    Variable            X 
-------------------------
          y2          100 
          z2          100 
          x2          0.5 
          x4          0.5 

  Constraint        Slack 
-------------------------
       A的需求量          100