I am modelling the logarithm of wages for 154 employees. The predictors are a binary matrix, where the rows are employees, the columns are work competencies, and the intersection is either 1 or 0 (i.e., whether the employee has the competency or not). I need to get integer coefficients for all 63 predictors. For the linear programming problem, I was given the objective function and the constraints: objective function: sum(n1+n2) --> min constraints: main constraint matrix: A x+E n1-E n2-y lambda - b = 0, where: A - binary matrix (154x63) x - integer coefficients which I need to take after modeling n1 - vector of excess logarithm of wages (positive difference between the predicted wage A*x and the real one) n2 - vector of logarithm deficiency (positive difference between the predicted wage A*x and the real one) y - logarithm of wages (154x1) lambda - some parameter, honestly I didn't understand from where I can take it, maybe from model, but I'm not sure (i think it is exogenous) (i took lambda.min from cv.glmnet() for the first time) b - similar to lambda E - maybe it is exogenous matrix x need to be integer x >= 1 I tried to write R code, but I don't know how linear programming works: obj <- c(1, 1) # Minimize sum(n1 + n2) A <- as.matrix(X) E <- matrix(rep(1), nrow = 154, ncol = 63) # just randomly for the first time b <- as.vector(offset_test) # just randomly for the first time y <- as.vector(log(grades_and_salaries$salary)) lambda <- as.vector(lambda.min) # just randomly mat <- cbind(A, E, -E, -y, rep(-lambda, nrow(A))) dir <- rep("==", nrow(A)) x.dir <- rep(">=", length(obj)) # x >= 1 mat <- rbind(mat, diag(length(obj))) dir <- c(dir, x.dir) rhs <- c(b, rep(1, length(obj))) lp.model <- lp("min", obj, mat, dir, rhs, all.int = TRUE) I got something strange, but I expected to get vector of integer coefficients for x. How can I solve my task in R with all constraints and get integer coefficients?

Reputation: 11

Integer Linear Programming in R (minimize the sum of differences between target variable and modeling results)

I am modelling the logarithm of wages for 154 employees. The predictors are a binary matrix, where the rows are employees, the columns are work competencies, and the intersection is either 1 or 0 (i.e., whether the employee has the competency or not).

I need to get integer coefficients for all 63 predictors.

For the linear programming problem, I was given the objective function and the constraints:

objective function: sum(n1+n2) --> min

constraints:

main constraint matrix: Ax+En1-En2-ylambda - b = 0, where:
- A - binary matrix (154x63)
- x - integer coefficients which I need to take after modeling
- n1 - vector of excess logarithm of wages (positive difference between the predicted wage A*x and the real one)
- n2 - vector of logarithm deficiency (positive difference between the predicted wage A*x and the real one)
- y - logarithm of wages (154x1)
- lambda - some parameter, honestly I didn't understand from where I can take it, maybe from model, but I'm not sure (i think it is exogenous) (i took lambda.min from cv.glmnet() for the first time)
- b - similar to lambda
- E - maybe it is exogenous matrix
x need to be integer
x >= 1

I tried to write R code, but I don't know how linear programming works:

obj <- c(1, 1)  # Minimize sum(n1 + n2)

A <- as.matrix(X)  
E <- matrix(rep(1), nrow = 154, ncol = 63)  # just randomly for the first time

b <- as.vector(offset_test)  # just randomly for the first time
y <- as.vector(log(grades_and_salaries$salary))  
lambda <- as.vector(lambda.min) # just randomly

mat <- cbind(A, E, -E, -y, rep(-lambda, nrow(A)))
dir <- rep("==", nrow(A))  

x.dir <- rep(">=", length(obj))  # x >= 1

mat <- rbind(mat, diag(length(obj)))
dir <- c(dir, x.dir)
rhs <- c(b, rep(1, length(obj)))

lp.model <- lp("min", obj, mat, dir, rhs, all.int = TRUE)

I got something strange, but I expected to get vector of integer coefficients for x.

How can I solve my task in R with all constraints and get integer coefficients?

Upvotes: 0

Answers (1)

Reinderien

Reputation: 15328

Running on a long list of assumptions and placeholder parameters,

import pandas as pd
import pulp
import numpy as np


rand = np.random.default_rng(seed=0)
m = 154
n = 63

A = rand.integers(low=0, high=2, size=(m, n))
x = pulp.LpVariable.matrix(
    name='x', indices=np.arange(n), cat=pulp.LpInteger,
    lowBound=1,
)  # (n, 1)
E = np.identity(m)  # (m, m)
n1 = pulp.LpVariable.matrix(
    name='n1', indices=np.arange(m), cat=pulp.LpContinuous,
    lowBound=0,
)  # (m, 1)
n2 = pulp.LpVariable.matrix(
    name='n2', indices=np.arange(m), cat=pulp.LpContinuous,
    lowBound=0,
)  # (m, 1)
salary = rand.uniform(low=100_000, high=300_000, size=m)
y = np.log(salary)  # (m, 1)
lambd = 1
b = offset_test = 20

prob = pulp.LpProblem('log_wage_model')
prob.setObjective(pulp.lpSum(n1) + pulp.lpSum(n2))

for i, lhs in enumerate(A@x + E@n1 - E@n2 - y*lambd - b):
    prob.addConstraint(
        name=f'constraint_{i}',
        constraint=lhs == 0,
    )

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal

df = pd.DataFrame({
    'A@x': A@x,
    'n1': n1,
    'n2': n2,
    'salary': salary,
    'y': y,
}).map(pulp.value)

pd.options.display.max_rows = m
print('Solution rows:')
print(df)

print('x =')
print([xi.value() for xi in x])

log_wage_model:
MINIMIZE
1*n1_0 + 1*n1_1 + 1*n1_10 + 1*n1_100 + 1*n1_101 + 1*n1_102 + 1*n1_103 + 1*n1_104 + 1*n1_105 + 1*n1_106 + 1*n1_107 + 1*n1_108 + 1*n1_109 + 1*n1_11 + 1*n1_110 + 1*n1_111 + 1*n1_112 + 1*n1_113 + 1*n1_114 + 1*n1_115 + 1*n1_116 + 1*n1_117 + 1*n1_118 + 1*n1_119 + 1*n1_12 + 1*n1_120 + 1*n1_121 + 1*n1_122 + 1*n1_123 + 1*n1_124 + 1*n1_125 + 1*n1_126 + 1*n1_127 + 1*n1_128 + 1*n1_129 + 1*n1_13 + 1*n1_130 + 1*n1_131 + 1*n1_132 + 1*n1_133 + 1*n1_134 + 1*n1_135 + 1*n1_136 + 1*n1_137 + 1*n1_138 + 1*n1_139 + 1*n1_14 + 1*n1_140 + 1*n1_141 + 1*n1_142 + 1*n1_143 + 1*n1_144 + 1*n1_145 + 1*n1_146 + 1*n1_147 + 1*n1_148 + 1*n1_149 + 1*n1_15 + 1*n1_150 + 1*n1_151 + 1*n1_152 + 1*n1_153 + 1*n1_16 + 1*n1_17 + 1*n1_18 + 1*n1_19 + 1*n1_2 + 1*n1_20 + 1*n1_21 + 1*n1_22 + 1*n1_23 + 1*n1_24 + 1*n1_25 + 1*n1_26 + 1*n1_27 + 1*n1_28 + 1*n1_29 + 1*n1_3 + 1*n1_30 + 1*n1_31 + 1*n1_32 + 1*n1_33 + 1*n1_34 + 1*n1_35 + 1*n1_36 + 1*n1_37 + 1*n1_38 + 1*n1_39 + 1*n1_4 + 1*n1_40 + 1*n1_41 + 1*n1_42 + 1*n1_43 + 1*n1_44 + 1*n1_45 + 1*n1_46 + 1*n1_47 + 1*n1_48 + 1*n1_49 + 1*n1_5 + 1*n1_50 + 1*n1_51 + 1*n1_52 + 1*n1_53 + 1*n1_54 + 1*n1_55 + 1*n1_56 + 1*n1_57 + 1*n1_58 + 1*n1_59 + 1*n1_6 + 1*n1_60 + 1*n1_61 + 1*n1_62 + 1*n1_63 + 1*n1_64 + 1*n1_65 + 1*n1_66 + 1*n1_67 + 1*n1_68 + 1*n1_69 + 1*n1_7 + 1*n1_70 + 1*n1_71 + 1*n1_72 + 1*n1_73 + 1*n1_74 + 1*n1_75 + 1*n1_76 + 1*n1_77 + 1*n1_78 + 1*n1_79 + 1*n1_8 + 1*n1_80 + 1*n1_81 + 1*n1_82 + 1*n1_83 + 1*n1_84 + 1*n1_85 + 1*n1_86 + 1*n1_87 + 1*n1_88 + 1*n1_89 + 1*n1_9 + 1*n1_90 + 1*n1_91 + 1*n1_92 + 1*n1_93 + 1*n1_94 + 1*n1_95 + 1*n1_96 + 1*n1_97 + 1*n1_98 + 1*n1_99 + 1*n2_0 + 1*n2_1 + 1*n2_10 + 1*n2_100 + 1*n2_101 + 1*n2_102 + 1*n2_103 + 1*n2_104 + 1*n2_105 + 1*n2_106 + 1*n2_107 + 1*n2_108 + 1*n2_109 + 1*n2_11 + 1*n2_110 + 1*n2_111 + 1*n2_112 + 1*n2_113 + 1*n2_114 + 1*n2_115 + 1*n2_116 + 1*n2_117 + 1*n2_118 + 1*n2_119 + 1*n2_12 + 1*n2_120 + 1*n2_121 + 1*n2_122 + 1*n2_123 + 1*n2_124 + 1*n2_125 + 1*n2_126 + 1*n2_127 + 1*n2_128 + 1*n2_129 + 1*n2_13 + 1*n2_130 + 1*n2_131 + 1*n2_132 + 1*n2_133 + 1*n2_134 + 1*n2_135 + 1*n2_136 + 1*n2_137 + 1*n2_138 + 1*n2_139 + 1*n2_14 + 1*n2_140 + 1*n2_141 + 1*n2_142 + 1*n2_143 + 1*n2_144 + 1*n2_145 + 1*n2_146 + 1*n2_147 + 1*n2_148 + 1*n2_149 + 1*n2_15 + 1*n2_150 + 1*n2_151 + 1*n2_152 + 1*n2_153 + 1*n2_16 + 1*n2_17 + 1*n2_18 + 1*n2_19 + 1*n2_2 + 1*n2_20 + 1*n2_21 + 1*n2_22 + 1*n2_23 + 1*n2_24 + 1*n2_25 + 1*n2_26 + 1*n2_27 + 1*n2_28 + 1*n2_29 + 1*n2_3 + 1*n2_30 + 1*n2_31 + 1*n2_32 + 1*n2_33 + 1*n2_34 + 1*n2_35 + 1*n2_36 + 1*n2_37 + 1*n2_38 + 1*n2_39 + 1*n2_4 + 1*n2_40 + 1*n2_41 + 1*n2_42 + 1*n2_43 + 1*n2_44 + 1*n2_45 + 1*n2_46 + 1*n2_47 + 1*n2_48 + 1*n2_49 + 1*n2_5 + 1*n2_50 + 1*n2_51 + 1*n2_52 + 1*n2_53 + 1*n2_54 + 1*n2_55 + 1*n2_56 + 1*n2_57 + 1*n2_58 + 1*n2_59 + 1*n2_6 + 1*n2_60 + 1*n2_61 + 1*n2_62 + 1*n2_63 + 1*n2_64 + 1*n2_65 + 1*n2_66 + 1*n2_67 + 1*n2_68 + 1*n2_69 + 1*n2_7 + 1*n2_70 + 1*n2_71 + 1*n2_72 + 1*n2_73 + 1*n2_74 + 1*n2_75 + 1*n2_76 + 1*n2_77 + 1*n2_78 + 1*n2_79 + 1*n2_8 + 1*n2_80 + 1*n2_81 + 1*n2_82 + 1*n2_83 + 1*n2_84 + 1*n2_85 + 1*n2_86 + 1*n2_87 + 1*n2_88 + 1*n2_89 + 1*n2_9 + 1*n2_90 + 1*n2_91 + 1*n2_92 + 1*n2_93 + 1*n2_94 + 1*n2_95 + 1*n2_96 + 1*n2_97 + 1*n2_98 + 1*n2_99 + 0
SUBJECT TO
constraint_0: n1_0 - n2_0 + x_0 + x_1 + x_10 + x_11 + x_12 + x_13 + x_14
 + x_15 + x_16 + x_17 + x_18 + x_19 + x_2 + x_21 + x_22 + x_25 + x_26 + x_28
 + x_29 + x_30 + x_33 + x_35 + x_45 + x_46 + x_47 + x_49 + x_50 + x_53 + x_54
 + x_55 + x_57 + x_58 + x_59 + x_60 + x_61 + x_62 + x_9 = 32.3426223434
...


Result - Optimal solution found

Objective value:                474.91974913
Enumerated nodes:               0
Total iterations:               11
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

Solution rows:
      A@x         n1        n2         salary          y
0    38.0   0.000000  5.657378  229262.369129  12.342622
1    33.0   0.000000  1.050450  154747.529521  11.949550
2    35.0   0.000000  2.437057  285770.012422  12.562943
3    32.0   0.555584  0.000000  283674.850854  12.555584
4    34.0   0.000000  1.883355  182890.970369  12.116645
5    36.0   0.000000  3.563959  251712.122649  12.436041
6    40.0   0.000000  7.532326  259801.652211  12.467674
7    32.0   0.284543  0.000000  216326.286444  12.284543
8    32.0   0.211185  0.000000  201025.020215  12.211185
9    34.0   0.000000  1.895780  180632.479017  12.104220
10   34.0   0.000000  1.743287  210388.906366  12.256713
11   29.0   2.991654  0.000000  161402.104187  11.991654
12   35.0   0.000000  3.125621  143541.359510  11.874378
13   39.0   0.000000  6.708636  217806.899571  12.291364
14   32.0   0.000000  0.403295  108738.869969  11.596705
15   35.0   0.000000  2.785452  201702.194806  12.214548
16   33.0   0.000000  0.503001  267533.168507  12.496999
17   33.0   0.000000  1.430002  105873.287671  11.569998
18   27.0   5.495064  0.000000  267016.144404  12.495064
19   34.0   0.000000  1.414815  292197.448173  12.585185
20   30.0   1.688665  0.000000  119212.703460  11.688665
21   26.0   5.648508  0.000000  114520.336704  11.648508
22   35.0   0.000000  3.428236  106060.364293  11.571764
23   37.0   0.000000  4.715235  216374.222582  12.284765
24   29.0   3.106623  0.000000  181067.076894  12.106623
25   28.0   4.366097  0.000000  234707.990334  12.366097
26   26.0   5.825136  0.000000  136644.271579  11.825136
27   38.0   0.000000  5.612143  239871.171418  12.387857
28   39.0   0.000000  6.872818  184828.264931  12.127182
29   36.0   0.000000  3.422853  289858.060771  12.577147
30   37.0   0.000000  4.620211  237943.686629  12.379789
31   33.0   0.000000  0.712445  216978.804476  12.287555
32   28.0   4.237749  0.000000  206436.608219  12.237749
33   31.0   1.116685  0.000000  182898.220970  12.116685
34   29.0   2.711578  0.000000  121975.752339  11.711578
35   33.0   0.000000  0.529253  260601.218073  12.470747
36   37.0   0.000000  4.656043  229568.457005  12.343957
37   38.0   0.000000  6.143451  141004.825949  11.856549
38   35.0   0.000000  3.412180  107777.081859  11.587820
39   29.0   2.704954  0.000000  121170.495922  11.704954
40   40.0   0.000000  8.458001  102949.979789  11.541999
41   28.0   4.588149  0.000000  293064.793331  12.588149
42   28.0   4.290188  0.000000  217550.833317  12.290188
43   31.0   1.117903  0.000000  183121.116609  12.117903
44   41.0   0.000000  8.930394  174487.026246  12.069606
45   34.0   0.000000  1.515953  264090.505973  12.484047
46   32.0   0.170466  0.000000  193003.982220  12.170466
47   33.0   0.000000  0.435772  286137.584139  12.564228
48   32.0   0.348152  0.000000  230533.600009  12.348152
49   34.0   0.000000  2.130800  142799.942126  11.869200
50   29.0   3.530806  0.000000  276732.277534  12.530806
51   28.0   4.420649  0.000000  247867.221786  12.420648
52   29.0   2.860175  0.000000  141516.958482  11.860175
53   31.0   1.312630  0.000000  222488.350025  12.312630
54   36.0   0.000000  3.768053  205242.470788  12.231947
55   36.0   0.000000  4.192889  134203.310824  11.807111
56   34.0   0.000000  1.784131  201968.932779  12.215869
57   25.0   7.535059  0.000000  277911.821637  12.535059
58   33.0   0.000000  0.517574  263662.781035  12.482426
59   20.0  12.529270  0.000000  276307.529033  12.529270
60   24.0   8.596200  0.000000  295433.861909  12.596200
61   32.0   0.609968  0.000000  299529.357840  12.609968
62   31.0   1.435923  0.000000  251682.269557  12.435923
63   26.0   6.599394  0.000000  296378.778383  12.599394
64   27.0   4.759453  0.000000  127957.402323  11.759453
65   25.0   7.154548  0.000000  189956.041211  12.154548
66   31.0   0.739899  0.000000  125479.677889  11.739899
67   34.0   0.000000  1.851949  188725.897408  12.148051
68   29.0   2.607723  0.000000  109943.649877  11.607723
69   28.0   3.551812  0.000000  103965.276043  11.551812
70   26.0   6.209559  0.000000  200698.541801  12.209559
71   33.0   0.000000  0.883416  182879.742710  12.116584
72   32.0   0.452657  0.000000  255929.462230  12.452657
73   33.0   0.000000  1.152469  139738.959296  11.847531
74   31.0   0.745966  0.000000  126243.297582  11.745966
75   31.0   0.739157  0.000000  125386.644699  11.739157
76   31.0   1.352831  0.000000  231614.891180  12.352831
77   34.0   0.000000  2.181739  135708.016320  11.818261
78   34.0   0.000000  2.043448  155834.916203  11.956552
79   35.0   0.000000  2.625445  236701.396495  12.374555
80   29.0   3.193201  0.000000  197442.217332  12.193201
81   29.0   2.945328  0.000000  154095.475075  11.945328
82   37.0   0.000000  4.779799  202845.746939  12.220201
83   29.0   2.774968  0.000000  129958.126942  11.774968
84   33.0   0.000000  0.993478  163819.764330  12.006522
85   25.0   7.478763  0.000000  262698.566072  12.478763
86   35.0   0.000000  2.469327  276695.510238  12.530673
87   33.0   0.000000  0.422917  289839.695523  12.577083
88   26.0   6.068594  0.000000  174310.572694  12.068594
89   32.0   0.000000  0.058682  153478.841506  11.941318
90   36.0   0.000000  4.417523  107202.721563  11.582477
91   30.0   2.321097  0.000000  224380.183706  12.321097
92   38.0   0.000000  5.978384  166311.150704  12.021616
93   30.0   1.596947  0.000000  108765.255675  11.596947
94   36.0   0.000000  3.757704  207377.495516  12.242296
95   29.0   2.807464  0.000000  134250.649582  11.807464
96   29.0   2.788956  0.000000  131788.753134  11.788956
97   27.0   4.606063  0.000000  109761.222428  11.606063
98   29.0   3.549336  0.000000  281908.003156  12.549336
99   32.0   0.250926  0.000000  209174.990413  12.250926
100  34.0   0.000000  2.261386  125318.482481  11.738614
101  28.0   3.771038  0.000000  129448.421070  11.771038
102  33.0   0.000000  0.430775  287570.942573  12.569225
103  28.0   3.653770  0.000000  115124.563443  11.653770
104  32.0   0.317166  0.000000  223499.806337  12.317166
105  29.0   2.571004  0.000000  105979.850267  11.571004
106  37.0   0.000000  4.919856  176335.560069  12.080144
107  30.0   2.572370  0.000000  288476.950924  12.572370
108  27.0   5.063236  0.000000  173379.131488  12.063236
109  28.0   4.149042  0.000000  188913.068814  12.149042
110  39.0   0.000000  7.354497  114176.760050  11.645503
111  30.0   1.646951  0.000000  114342.260909  11.646952
112  28.0   4.277224  0.000000  214748.661155  12.277224
113  36.0   0.000000  3.419639  290791.265371  12.580361
114  31.0   1.518807  0.000000  273431.609899  12.518807
115  32.0   0.346233  0.000000  230091.695188  12.346233
116  28.0   4.597144  0.000000  295712.813714  12.597144
117  35.0   0.000000  2.864618  186349.952387  12.135382
118  32.0   0.582983  0.000000  291554.622472  12.582983
119  30.0   2.229906  0.000000  204823.881443  12.229906
120  32.0   0.000000  0.367599  112690.586456  11.632401
121  34.0   0.000000  1.830026  192908.986342  12.169974
122  34.0   0.000000  1.410835  293362.630228  12.589165
123  27.0   5.119152  0.000000  183349.874215  12.119151
124  37.0   0.000000  4.398223  297086.063654  12.601777
125  40.0   0.000000  7.911671  177784.750619  12.088329
126  37.0   0.000000  5.206957  132328.502235  11.793043
127  33.0   0.000000  1.075586  150906.229495  11.924414
128  29.0   2.904558  0.000000  147939.422941  11.904558
129  33.0   0.000000  0.714688  216492.630817  12.285312
130  34.0   0.000000  2.210345  131880.932157  11.789655
131  35.0   0.000000  2.946624  171677.993782  12.053376
132  26.0   6.591676  0.000000  294100.319926  12.591676
133  32.0   0.000000  0.279193  123106.697198  11.720807
134  32.0   0.348017  0.000000  230502.475417  12.348017
135  37.0   0.000000  5.113488  145293.624324  11.886512
136  33.0   0.000000  0.893438  181056.077762  12.106562
137  31.0   0.657601  0.000000  115566.445335  11.657601
138  33.0   0.000000  0.567016  250943.712535  12.432984
139  31.0   1.505715  0.000000  269875.220048  12.505715
140  31.0   1.435954  0.000000  251690.151225  12.435954
141  35.0   0.000000  2.981896  165728.138731  12.018104
142  29.0   3.246344  0.000000  208218.626320  12.246344
143  36.0   0.000000  4.000142  162731.610654  11.999858
144  33.0   0.000000  0.402298  295877.908802  12.597702
145  34.0   0.000000  1.869567  185430.084963  12.130433
146  35.0   0.000000  2.517524  263675.966099  12.482476
147  34.0   0.000000  2.321838  117967.243305  11.678162
148  31.0   1.575606  0.000000  289411.867920  12.575606
149  34.0   0.000000  1.468080  277040.844396  12.531920
150  36.0   0.000000  3.619432  238129.090308  12.380568
151  38.0   0.000000  5.812084  196401.410072  12.187916
152  29.0   3.371170  0.000000  235901.626528  12.371170
153  41.0   0.000000  8.504428  267151.750313  12.495572
x =
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

Or if you make lambda and b degrees of freedom:

lambd = pulp.LpVariable(name='lambda', cat=pulp.LpContinuous, lowBound=0.5)
b = pulp.LpVariable(name='b', cat=pulp.LpContinuous)
# ...
print('lambda =', lambd.value())
print('b =', b.value())

Objective value:                474.13786593
Enumerated nodes:               4
Total iterations:               965
Time (CPU seconds):             0.15
Time (Wallclock seconds):       0.16

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.15   (Wallclock seconds):       0.16

lambda = 1.1403868
b = 18.633747

Upvotes: 0

Integer Linear Programming in R (minimize the sum of differences between target variable and modeling results)

Answers (1)

Related Questions