The problem being solved is finding the truss with the least weight, exactly done as on this website: https://www.layopt.com/truss/. This method is also called the ground structure method. I am aiming to add some functionallity to the method through an MILP optimization. Initially, I started with the cvxpy
Linear Programming (LP) solver but since it can't solve MILP I am now using the IBM ILOG CPLEX solver with the python module docplex. Unfortunately I am bumping into an problem, but first in short how the optimization works is:
- Initialize: Make nodes, place a load (
f
) and supports (dof
) - Make the potential members of the truss, with full connectivity all the nodes will be connected with each other. These potential members will form the truss and every member has two variables
a
(area) andq
(internal force) and a lengthl
. - Calculate the nodal equilibrium matrix
B
, which is guaranteed sparce and symmetric. - Build the LP-model, which is done via calling the
Model()
object and start solving it.
The formulation is: \begin{array}{ll} \text{minimize} & V = l^T a \\ \text{subject to}& B q = f \\ & q \le \sigma a \\ & q \ge -\sigma a \\ & a \ge 0 \\ \end{array} In code:
"build the model"
model = Model()
a = model.continuous_var_list(len(members), name='a', lb=0, ub=1)
q = []
Bcons = []
for k, fk in enumerate(f):
qk = model.continuous_var_list(len(members), name=f'q{k}', lb=-10 ** 10, ub=10 ** 10)
Bconsk = model.add_constraints(
sum(B[i, j] * qk[j] for j in range(len(qk))) == fk[i] for i in range(len(dof)) if dof[i] != 0)
model.add_constraints(qk[i] <= sigma * a[i] for i in range(len(qk)))
model.add_constraints(qk[i] >= -sigma * a[i] for i in range(len(qk)))
q.append(qk)
Bcons.append(Bconsk)
model.set_objective('min', sum(a[i] * l[i] for i in range(len(a))))
model.print_information()
"solve model and update results"
sol = model.solve()
print(model.solve_details)
For small problem sizes (less than 1000 members) the LP converges in the CPLEX solver and works perfectly. However for larger problem sizes (e.g. more than 10 000 members) the problem is 'optimal with unscaled infeasibilities'. What does this mean and how can I solve this?
What I tried so far:
- Solving the same problem in cvxpy. This module can solve up to 120 000 members without any problems.. So I actually would expect that the professional solver can also do this.
- Found Coping with an Ill-Conditioned Problem or Handling Unscaled Infeasibilities, but I cant find any scaling parameters in the module
docplex
.
..
Results:
Nodes: 231 Members: 16290
Model: docplex_model1
- number of variables: 32580
- binary=0, integer=0, continuous=32580
- number of constraints: 33038
- linear=33038
- parameters: defaults
- objective: minimize
- problem type is: LP
status = optimal with unscaled infeasibilities
time = 131.297 s.
problem = LP