У меня есть эта функция без штрафов в целых числах решённая. Но я не пойму, что ещё нужно дописать в функцию, чтобы выполнялись штрафы и чтобы это было линейным также.
Замечательное свойство транспортной задачи состоит в том, что матрица там абсолютно унимодулярна (т.е. определитель равен

), и поэтому решение при целочисленных объёмах поставок тоже будет целочисленным. Введение дополнительных ограничений это нарушает, даже линейных. А тут вообще нелинейные. Притом дискретные (в классической ТЗ нам дико везёт, дискретность получается "сама собой", вообще же говоря это чудовищное усложнение).
Кстати, штрафы за перевозку по определённому маршруту (

) или за привлечение определённого поставщика? (далее принимаю первое).
Можно рассматривать задачу, как задачу ЦЛП общего вида, добавив к nm целочисленным переменным

, соответствующим перевозкам от i-того поставщика к j-тому потребителю столько же "штрафных"

, так же целочисленных, связанных с основными неравенствами

, где M - достаточно большое число (т.е. игрек равен нулю, только если х равен нулю), в ЦФ игреки входят с коэффициентами, равными штрафам. Размерность задачи огромная, и целочисленное решение может потребовать много времени.
У Триуса "Задачи математического программирования транспортного типа" рассмотрен вариант нелинейных функций потерь, но в данной задаче это не более чем аппроксимация. Может, смотреть в сторону сетевых алгоритмов. Но в любом случае такая функция, с "порогом", сильно усложняет дело.