Why the result is infeasible?

Suppose I have a set of jobs, each needing one mold; changing from mold A to mold B has a cost. I want to search for the job sequence that minimizes the mold change cost.

job mold

Mold shift cost

My code is as follows:

%pip install -q amplpy

SOLVER = "highs"

from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

import pandas as pd

jobmold = pd.DataFrame(
    {
        "Job1": {"mold": 1},
        "Job2": {"mold": 2},
        "Job3": {"mold": 1},
        "Job4": {"mold": 3},
        "Job5": {"mold": 2},
        "Job6": {"mold": 3},
        "Job7": {"mold": 2},
    }
).T

moldshiftcost = pd.DataFrame(
    {
        "mold1": {"mold1": 0, "mold2": 2, "mold3": 4},
        "mold2": {"mold1": 2, "mold2": 0, "mold3": 4},
        "mold3": {"mold1": 7, "mold2": 1, "mold3": 0}
    }
).T

display(jobmold)
display(moldshiftcost)

%%writefile LPanelProduction.mod
set JOBS;
set MOLDS;
set jobpairs within {JOBS,JOBS};
set moldpairs within {MOLDS,MOLDS};
set SHIFTJOBS;

param jobshiftcost{jobpairs};
param moldshiftcost{moldpairs};
param jn;
param mn;
param shiftnum;

var sequence{JOBS} integer >=0,<jn

minimize obj:sum{k in SHIFTJOBS, i in JOBS, j in JOBS}(if sequence[k]==i and sequence[k+1]==j then 1 else 0)*jobshiftcost[i,j];

subj to no_overlap{i in JOBS, j in JOBS}:(sequence[i]> sequence[j] or sequence[i]<sequence[j]);

def Lpanelproduction(jobmold, moldshiftcost):
    model = AMPL()
    model.read("LPanelProduction.mod")

    jn = len(jobmold)
    JOBS = range(jn)
    SHIFTJOBS = range(jn - 1)
    shiftnum = len(SHIFTJOBS)
    jobpairs = {(i, j) for i in JOBS for j in JOBS}

    model.param["jn"] = jn
    model.set["JOBS"] = JOBS
    model.set["SHIFTJOBS"] = SHIFTJOBS
    model.set["jobpairs"] = jobpairs
    model.param["shiftnum"] = shiftnum

    mn = len(moldshiftcost)
    MOLDS = range(mn)
    moldpairs = {(i, j) for i in MOLDS for j in MOLDS}

    model.param["mn"] = mn
    model.set["MOLDS"] = MOLDS
    model.set["moldpairs"] = moldpairs

    moldshiftcost = moldshiftcost.values
    model.param["moldshiftcost"] = moldshiftcost

    jobmold = jobmold.values

    jobshiftcost = {(i, j): moldshiftcost[(k, l)].tolist() for (i, j) in jobpairs for (k, l) in moldpairs if
                    (jobmold[i] - 1) == k and (jobmold[j] - 1) == l}
    model.param["jobshiftcost"] = jobshiftcost

    return model

y = Lpanelproduction(jobmold, moldshiftcost)
y.solve(solver=SOLVER)
assert y.solve_result == "solved", y.solve_result

y.display("sequence")

Why can't I obtain a feasible solution?

Your issue is that you have a syntax error when declaring the variables, you forgot the semicolon at the end.

Besides that, I suggest replacing the no_overlap constraint with this cleaner alldiff:

subject to visit_all_jobs:
    alldiff{j in JOBS} sequence[j]; # since every job in the sequence is different, this way you look for every permutation

However, when you do this, HiGHS solver can't figure out bounds for your variables because you used a strict bound (seq[k] < jn). In this case, since sequence is an integer, you have to use (seq[k] <= jn-1).

The objective is determined by the costs from moving between molds, that you can write cleaner:

minimize cost:
    sum{k in SHIFTJOBS, i in JOBS, j in JOBS}
        (if sequence[k]==i and sequence[k+1]==j then jobshiftcost[i,j]);

If Ampl supported variables in indexes, you could do something easier like this:

minimize cost:
    sum{k in SHIFTJOBS} jobshiftcost[seq[k],seq[k+1]];

But what you have is quite compact already.

Final model:

set JOBS;
set MOLDS;
set jobpairs within {JOBS,JOBS};
set moldpairs within {MOLDS,MOLDS};
set SHIFTJOBS;

param jobshiftcost{jobpairs};
param moldshiftcost{moldpairs};
param jn;
param mn;
param shiftnum;

var sequence{JOBS} integer >=0,<=jn-1;

minimize cost:
    sum{k in SHIFTJOBS, i in JOBS, j in JOBS}
        if sequence[k]==i and sequence[k+1]==j then jobshiftcost[i,j];

subject to visit_all_jobs:
    alldiff{k in JOBS} sequence[k]; # since every sequence[k] is different to sequence[k'], this way you look for every permutation

The solution I get for your example is:

y = Lpanelproduction(jobmold, moldshiftcost)
y.solve(solver=SOLVER, verbose=True)
assert y.solve_result == "solved", y.solve_result
sol = y.var['sequence'].to_dict()
print(sol)
> {0: 3, 1: 5, 2: 1, 3: 4, 4: 6, 5: 2, 6: 0}