How to solve multiple multivariate equation systems with constraints

Question

I am trying to solve a blending problem with a system of 3 equations and I have 3 objectives to reach for, or try to get the values as close as posible for the three of them: The equations are:

def sat (c,s,a,f):
  return (100*c)/(2.8*s+1.18*a+0.65*f) #For this I need sat = 98.5

def ms (s,a,f):
  return s/(a+f)            #For this I need ms = 2.5
 
def ma (a,f):
  return (a/f)              #For this I need ms = 1.3


#The total mix ratio:
r1+r2+r3+r4+r5+r6 = 1

material_1:
c = 51.29
s = 4.16
a = 0.97
f = 0.38

material_2:
c = 51.42
s = 4.16
a = 0.95
f = 0.37

material_3:

c = 6.88
s = 63.36
a = 13.58
f = 3.06

material_4:
c = 32.05
s = 1.94
a = 0.0
f = 0.0

material_5:
c = 4.56
s = 21.43
a = 3.82
f = 52.28

material_6:
c = 0.19
s = 7.45
a = 4.58
f = 0.42

#The aproximate values I am trying to find are around: 

0.300 <= r1 <= 0.370 
0.300 <= r2 <= 0.370 
0.070 <= r3 <= 0.130 
0.005 <= r4 <= 0.015 
0.010 <= r5 <= 0.030 
0.110 <= r6 <= 0.130 

So how can I calculate the value for every ratio "r" in order to get the closets values to the objectives for the 3 equations? I looked on some optimizers but as I am new with them I still can not understand how to set up the problem, the equations and constraints into them.

Answer 1

Ok this is what i have done and works fine:

#I call the values from a pandas DF:
    c1 = df.at[0, 'MAX']
    c2 = df.at[4, 'MAX']
    c3 = df.at[8, 'MAX']
    c5 = df.at[12, 'MAX']
    c6 = df.at[16, 'MAX']
    c7 = df.at[20, 'MAX']
    s1 = df.at[1, 'MAX']
    s2 = df.at[5, 'MAX']
    s3 = df.at[9, 'MAX']
    s5 = df.at[13, 'MAX']
    s6 = df.at[17, 'MAX']
    s7 = df.at[21, 'MAX']
    a1 = df.at[2, 'MAX']
    a2 = df.at[6, 'MAX']
    a3 = df.at[10, 'MAX']
    a5 = df.at[14, 'MAX']
    a6 = df.at[18, 'MAX']
    a7 = df.at[22, 'MAX']
    f1 = df.at[3, 'MAX']
    f2 = df.at[7, 'MAX']
    f3 = df.at[11, 'MAX']
    f5 = df.at[15, 'MAX']
    f6 = df.at[19, 'MAX']
    f7 = df.at[23, 'MAX']


r1 = cp.Variable()
r2 = cp.Variable()
r3 = cp.Variable()
r5 = cp.Variable()
r6 = cp.Variable()
r7 = 12.5

#Objectives
ma = 1.3
ms = 2.50
lsf = 98.5


delta1 =(ms*((r1*a1+r2*a2+r3*a3+r5*a5+r6*a6+r7*a7)+(f1*r1+f2*r2+f3*r3+f5*r5+f6*r6+f7*r7))-(r1*s1+r2*s2+r3*s3+r5*s5+r6*s6+r7*s7))**2
delta2 =(ma*(f1*r1+f2*r2+f3*r3+f5*r5+f6*r6+f7*r7)-(r1*a1+r2*a2+r3*a3+r5*a5+r6*a6+r7*a7))**2
delta3 =((lsf*(2.8*(r1*s1+r2*s2+r3*s3+r5*s5+r6*s6+r7*s7)+1.18*(r1*a1+r2*a2+r3*a3+r5*a5+r6*a6+r7*a7)+0.65*(f1*r1+f2*r2+f3*r3+f5*r5+f6*r6+f7*r7))-100*(r1*c1+r2*c2+r3*c3+r5*c5+r6*c6+r7*c7)))**2

objective = cp.Minimize(delta1+delta2+delta3)

constraints = [r1-r2 == 0, #I added this to make r1=r2.
               r1>= 0.20,
               r3>=0,      #I added these to make it non negative.
               r5>=0,
               r5<=0.008,
               r6>=0,
               r1+r2+r3+r5+r6+r7 == 1]
problem = cp.Problem(objective,constraints)
problem.solve()
print(r1.value,r2.value,r3.value,r5.value,r6.value)
print(problem.status)

Once again i want to thank you for your help guys.

Maybe you know how I can improve the code for get the variable values, maybe there is and example of using a for loop to get the values instead of put it directly from the DF for each one, the DF looks like this:

    DATO    MAX
0   c1  51.95000    
1   s1  3.07000 
2   a1  0.83000 
3   f1  0.31000 
4   c2  52.26000    
5   s2  2.82000 
6   a2  0.75000
... 

Answer 2

So this is a simple/trivial example to show the intent that was mentioned in comment to minimize the square of errors... Instead of using a constraint to pin a value to an exact outcome, we let the solver find the best outcome that minimizes the square of the error where error = value - target. I think what I've written below is fairly clear. CVXPY likes to work in the linear algebra realm, and I'm sure this could be converted into vector / matrix format, but the concept is to remove constraints and let the solver figure out the best combo. Obviously, if there are hard constraints, those need to be added, but note that I've just made an example with 2 of your 3 targets (with some trivial math) and moved it into the objective.

Your problem with 3 constraints that aren't simultaneously satisfiable is probably a candidate for a conversion like this...

import cvxpy as cp

r1 = cp.Variable()
r2 = cp.Variable()

ma = 2.5
ms = 3.4

delta_1 = (r1 + r2 - ma)**2      # diff from r1 + r2 and ma
delta_2 = (3*r1 + 2*r2 - ms)**2  # diff from 3r1 + 2r2 and ms

prob = cp.Problem(cp.Minimize(delta_1 + delta_2))

prob.solve()
print(prob.value)
print(r1.value, r2.value)

Output

9.860761315262648e-31
-1.6000000000000014 4.100000000000002

Answer 3

I guess I made it, of course the code is awful but I will try to make it look better later. I added the cost of the components so I can give a function to "minimize", of course this is becouse I know the aproximated material ratio so it guide the solver to it.

I will post the code for it:

c1 = 51.42
c2 = 51.42
c3 = 6.88
c5 = 32.05
c6 = 4.56
c7 = 0.19

s1 = 4.16
s2 = 4.16
s3 = 63.36
s5 = 1.94
s6 = 21.43
s7 = 7.45

a1 = 0.97
a2 = 0.95
a3 = 13.58
a5 = 0.0
a6 = 3.82
a7 = 4.58

f1 = 0.38
f2 = 0.37
f3 = 3.06
f5 = 0.0
f6 = 52.28
f7 = 0.42

r7 = 0.125

r1 = cp.Variable()
r2 = cp.Variable()
r3 = cp.Variable()
r5 = cp.Variable()
r6 = cp.Variable()

#Costos
caliza = 10
arcilla = 20
hierro = 170
yeso = 80


objective = cp.Minimize(r1*caliza+r2*caliza+r3*arcilla+r5*yeso+r6*hierro)
constraints = [
               r1-r2 == 0,
               r1>= 0.20,
               r1<= 0.40,
               r3<=0.14,
               r3>=0.06,
               r5>=0.001,
               r5<=0.008,
               r6>=0.01,
               r6<=0.03,
               2.5*((r1*a1+r2*a2+r3*a3+r5*a5+r6*a6+r7*a7)+(f1*r1+f2*r2+f3*r3+f5*r5+f6*r6+f7*r7))-(r1*s1+r2*s2+r3*s3+r5*s5+r6*s6+r7*s7)==0,
               (98.5*(2.8*(r1*s1+r2*s2+r3*s3+r5*s5+r6*s6+r7*s7)+1.18*(r1*a1+r2*a2+r3*a3+r5*a5+r6*a6+r7*a7)+0.65*(f1*r1+f2*r2+f3*r3+f5*r5+f6*r6+f7*r7))-100*(r1*c1+r2*c2+r3*c3+r5*c5+r6*c6+r7*c7)) == 0,
               #1.3*(f1*r1+f2*r2+f3*r3+f5*r5+f6*r6+f7*r7)-(r1*a1+r2*a2+r3*a3+r5*a5+r6*a6+r7*a7) == 0,
               r1+r2+r3+r5+r6+r7 == 1]
problem = cp.Problem(objective,constraints)
problem.solve()
print(r1.value,r2.value,r3.value,r5.value,r6.value)
print(problem.status)

This gives me the result:

0.3644382497863931 0.3644382497863931 0.12287226775076901 0.0009999999955268117 0.022251232680917873
optimal

Anyways, the only way to make a feasible result is to only consider 2 of the three constraints functions, becouse the components cant reach the 3 of them and this indicates that I need to check de material components before I try to reach the 3 constraints (wich were the sat, ma and ms). Now I will try to make the code better using pandas so i can get the material components with somekind of for loop and laso use it for the ratios. Thank you so much for your help👍.