ILP runtime seems to be linear? I have a variation the shortest path problem, formulated as an ILP. 
The system model is as follows:


There is a connected digraph consisting of 20 nodes, with each link having an associated weight delay which is a float. The graph is modeled as a python networkx object as follows:

g = nx.DiGraph(nx.barabasi_albert_graph(20, 2))
for u, v in g.edges:
    g[u][v]['delay'] = random.uniform(50, 1000)  # uS

There are multiple source-destination pairs and each pair has an associated 'demand' in terms of delay. The problem is modeled in python PuLP and solved using Gurobi solver.

# get cost
delay_d = nx.get_edge_attributes(self.graph, 'delay')


# get links
links = []
for i, j in self.graph.edges:
    links.append((i, j))

# instantiate the problem
prob = pulp.LpProblem("Shortest Path Problem", pulp.LpMinimize)

# create binary variables to state a link is chosen on shortest path
var_dict = {}
for flow_name in self.list_of_flows:
    for (i, j) in links:
        var = pulp.LpVariable('x_(%s,%s)_%s' % (i, j, flow_name), cat=pulp.LpBinary)
        var_dict[(i, j, flow_name)] = var

# create slack variables for every flow
delay_slack = {}
for flow_name in self.list_of_flows:
    delay_slack[flow_name] = pulp.LpVariable('d_%s' % flow_name, lowBound=0.0, cat=pulp.LpContinuous)

# formulate the objective
prob += pulp.lpSum([pulp.lpSum([delay_d[(i, j)] * var_dict[(i, j, flow_name)] for (i, j) in links])
                    for flow_name in self.list_of_flows ]) \
        + pulp.lpSum([1000 * delay_slack[flow_name] for flow_name in self.list_of_flows])

# formulate the constraints
for flow_name, flow in self.list_of_flows.items():
    # conservation of flow constraints
    for node in self.graph.nodes:
        if node == flow['source']:
            prob += pulp.lpSum([var_dict[(i, j, flow_name)] for (i, j) in links if i == node]) - \
                    pulp.lpSum([var_dict[(j, i, flow_name)] for (j, i) in links if i == node]) == 1
        elif node == flow['target']:
            prob += pulp.lpSum([var_dict[(i, j, flow_name)] for (i, j) in links if i == node]) - \
                    pulp.lpSum([var_dict[(j, i, flow_name)] for (j, i) in links if i == node]) == -1
        else:
            prob += pulp.lpSum([var_dict[(i, j, flow_name)] for (i, j) in links if i == node]) - \
                    pulp.lpSum([var_dict[(j, i, flow_name)] for (j, i) in links if i == node]) == 0

    # delay constraint
    prob += pulp.lpSum([var_dict[(i, j, flow_name)] * delay_d[(i, j)]
                        for (i, j) in links]) - delay_slack[flow_name] &lt;= flow['qos']['delay']

# solve the optimization problem
prob.solve(pulp.GUROBI_CMD(msg=0))



There are no issues with the solution and correct results are obtained.
However, when I calculated the runtime, the results seem to be linear.



Now, I know that if an ILP satisfies the property of unimodularity, then the LP relaxation gives integer solution. Thus such an ILP can be solved in polynomial time using methods such as simplex. 

However, if we look at the delay constraints, the coefficients of LHS are not integers, nor is the RHS integer. Thus it seems unimodularity property is not satisfied. Hence, shouldn't the runtime be exponential?

Can anyone explain why the runtime graph is linear?
I will try to answer to your question with an example.
Consider the following time complexity: $f(n)=2^\frac{n}{1000}$. 

It is clearly exponential in $n$, but look at the plot.

For $n \in [0,80]$:



From this perspective it looks constant, right?

For $n \in [0, 1.5 \cdot 10^4]$:



From this perspective everything changes.

What I want to show, is that inferring the complexity just looking at the runtime with different inputs size is not possibile.