# Is there any way to model this situation in integer programming?

I have a binary variable $z_{ij}$. Some constraints are given by: $$a_j\geq b_i z_{ij},$$ where $a_j$ and $b_i$ are given nonnegative inputs BUT $a_j$ changes depending on the solution $z_{ij}$. I mean, initially (before solving the problem), I know the values of $a_j$ for all $j$ and $$\mathbf{a}=f(\mathbf{Z}),$$ where $\mathbf{a}=(a_1,\ldots,a_n)$ and $\mathbf{Z}=[z_{ij}]$. Also, the function $f(\cdot)$ is not analytically (explicitly) given, i.e., I have an algorithm that calculates $f(\cdot)$ but I don't have the exact formula for it.

Can I still model this problem using integer programming? If not, what tool can I use?

If $a$ is dependent on the decision variables, then it's not an input. You may be using $a$ to represent some quantity that exists both before and after the optimisation process (and is changed by optimisation) but in that case, the "before optimisation" and "after optimisation" versions are best treated as separate entities: one is an input, the other is dependent on decision variables.

If $a$ can be expressed as a linear function of $z$ then it's easy to convert to a linear problem, but from your question I gather this isn't the case. That means the general problem that you're describing is nonlinear, and potentially very hard if $f()$ is effectively a black-box function.

That said, sometimes it's possible to tackle such problems in an integer framework, depending on the nature of the problem and how squeamish you are :-)

Option 1: treat the initial values of $a$ as fixed inputs, base your constraints on these values, and solve. When you get a solution, update the values of $a$ used in the constraints according to that value, and re-solve until the solution converges to something where the values of $Z$ and the values of $a$ derived from those values satisfy the constraints given.

Option 2: Using whatever method seems appropriate (e.g. by plugging in those initial values of $a$), find an initial estimate of what the solutions for $z$ will be.

Now find the partial derivatives of $f(Z)$ in the neighbourhood of this estimate. Ideally you'd do this by differentiating $f$ and defining a new function $g(z) = f'(z)$ but if $f$ is a black box you'll have to do it the hard way, by calculating values for nearby points and then using those to estimate a gradient.

From this you can create a function $h(x)$ which is linear, and approximates $f$ in the neighbourhood of your original solution. You can then replace $a$ with $g(Z)$ in your constraints to get a linear constraint that approximates the nonlinear constraint. You can then use the same "solve, update estimates, repeat until converged" approach as for Option 1 above.

Option 3: lazy constraint generation: ignore constraints, solve, check whether solution satisfies constraints. If it doesn't, identify a combination of values that leads to constraint violation, add a constraint that forbids this particular combination, and re-solve.

Whether any of these will work for your particular problem depends on a lot of specifics. In theory, Option 3 should eventually give you an optimal solution (since your search space is finite) but run-time may be prohibitive. Options 1 and 2 are likely to struggle if $f$ isn't smooth; for some problems they may fail to converge to a legitimate solution, and for others they may deliver a valid but suboptimal solution.