""" City Resource Optimization -> CP-SAT (Google OR-Tools) Maximise Electrum * Brass * Steel after the gains of step 5. The objective is a product of three variables, so this is a constraint- programming / nonlinear problem, hence CP-SAT (cp_model) rather than the LP/MIP solver. Read the comments at: - PARAMETERS (initial pools + arrival schedule + bonus mode) - "OBJECTIVE IS SET HERE" (the product being maximised -- tweak freely) """ from ortools.sat.python import cp_model import printer # ====================================================================== # PARAMETERS -- edit these # ====================================================================== # Starting resource pools at the start of step 1: (Electrum, Brass, Steel, Capital) INITIAL = (3, 3, 3, 3) # Arrival schedule. Key = step (1..5), value = list of city types that arrive # at the START of that step. Types: 'H' Hub, 'F' Foundry, 'M' Metropolis, 'N' Monument. # Total cities across all steps must be <= 7. Arriving cities act that same step. ARRIVALS = { 1: ["H", "F", "H"], 2: ["F"], 3: [], 4: ["N"], 5: [], } # Collect Bonus (b) adds +1 to whatever a Collect gives. On a foundry that is # +1 of the resource collected (i.e. the chosen vat's value + 1). This is the # same uniform rule as Hub (+1 Capital) and Metropolis (+1 resource pick). NUM_STEPS = 5 MAX_RES = ( 200 # upper bound on any resource pool (raise if you expect more; affects speed) ) MAX_VAT = 15 # upper bound on a foundry vat value (1 + 2*nsteps is plenty) # ====================================================================== # MODEL # ====================================================================== def solve( initial=INITIAL, arrivals=ARRIVALS, max_res=MAX_RES, max_vat=MAX_VAT, time_limit=60.0, num_workers=8, verbose=True, ): # ---- build the city list ----------------------------------------- cities = [] # list of (arrival_step, arrival_type) for s in range(1, NUM_STEPS + 1): for typ in arrivals.get(s, []): cities.append((s, typ)) N = len(cities) assert N <= 7, f"At most 7 cities allowed, got {N}" m = cp_model.CpModel() def AND(a, b): """Boolean AND of two 0/1 vars, returned as a new 0/1 var.""" c = m.NewBoolVar("") m.AddMultiplicationEquality(c, [a, b]) return c # ---- state variables --------------------------------------------- # Indexed [i][t]; t in 1..NUM_STEPS+1 for state (t=NUM_STEPS+1 == "after step 5"). isH, isF, isM, isMon, present = {}, {}, {}, {}, {} hasA, hasB, hasD = {}, {}, {} vE, vB, vS = {}, {}, {} # foundry vats at start of step t # action variables, t in 1..NUM_STEPS col, ua, ub, ud = {}, {}, {}, {} # collect / upgrade a,b,d ow = {} # overwork (global limit: 1 per step) rH, rF, rM = {}, {}, {} # renovate -> Hub/Foundry/Metro cvE, cvB, cvS = {}, {}, {} # foundry: which vat collected (normal collect) owcvE, owcvB, owcvS = {}, {}, {} # foundry: which vat collected (overwork) mE, mB, mS, mC = {}, {}, {}, {} # metropolis: resources picked (+1 each) owmE, owmB, owmS, owmC = {}, {}, {}, {} # metropolis: resources picked (overwork) for i in range(N): a_step, a_type = cities[i] for t in range(1, NUM_STEPS + 2): isH[i, t] = m.NewBoolVar(f"isH_{i}_{t}") isF[i, t] = m.NewBoolVar(f"isF_{i}_{t}") isM[i, t] = m.NewBoolVar(f"isM_{i}_{t}") isMon[i, t] = m.NewBoolVar(f"isMon_{i}_{t}") present[i, t] = m.NewBoolVar(f"present_{i}_{t}") hasA[i, t] = m.NewBoolVar(f"hasA_{i}_{t}") hasB[i, t] = m.NewBoolVar(f"hasB_{i}_{t}") hasD[i, t] = m.NewBoolVar(f"hasD_{i}_{t}") vE[i, t] = m.NewIntVar(0, max_vat, f"vE_{i}_{t}") vB[i, t] = m.NewIntVar(0, max_vat, f"vB_{i}_{t}") vS[i, t] = m.NewIntVar(0, max_vat, f"vS_{i}_{t}") # presence: present from arrival step onward, persists m.Add(present[i, t] == (1 if t >= a_step else 0)) # exactly one type iff present m.Add(isH[i, t] + isF[i, t] + isM[i, t] + isMon[i, t] == present[i, t]) for t in range(1, NUM_STEPS + 1): col[i, t] = m.NewBoolVar(f"col_{i}_{t}") ua[i, t] = m.NewBoolVar(f"ua_{i}_{t}") ub[i, t] = m.NewBoolVar(f"ub_{i}_{t}") ud[i, t] = m.NewBoolVar(f"ud_{i}_{t}") ow[i, t] = m.NewBoolVar(f"ow_{i}_{t}") rH[i, t] = m.NewBoolVar(f"rH_{i}_{t}") rF[i, t] = m.NewBoolVar(f"rF_{i}_{t}") rM[i, t] = m.NewBoolVar(f"rM_{i}_{t}") cvE[i, t] = m.NewBoolVar(f"cvE_{i}_{t}") cvB[i, t] = m.NewBoolVar(f"cvB_{i}_{t}") cvS[i, t] = m.NewBoolVar(f"cvS_{i}_{t}") owcvE[i, t] = m.NewBoolVar(f"owcvE_{i}_{t}") owcvB[i, t] = m.NewBoolVar(f"owcvB_{i}_{t}") owcvS[i, t] = m.NewBoolVar(f"owcvS_{i}_{t}") mE[i, t] = m.NewIntVar(0, 3, f"mE_{i}_{t}") mB[i, t] = m.NewIntVar(0, 3, f"mB_{i}_{t}") mS[i, t] = m.NewIntVar(0, 3, f"mS_{i}_{t}") mC[i, t] = m.NewIntVar(0, 3, f"mC_{i}_{t}") owmE[i, t] = m.NewIntVar(0, 6, f"owmE_{i}_{t}") owmB[i, t] = m.NewIntVar(0, 6, f"owmB_{i}_{t}") owmS[i, t] = m.NewIntVar(0, 6, f"owmS_{i}_{t}") owmC[i, t] = m.NewIntVar(0, 6, f"owmC_{i}_{t}") # ---- per-step gain/cost accumulators (linear expressions) --------- gain_E = {t: [] for t in range(1, NUM_STEPS + 1)} gain_B = {t: [] for t in range(1, NUM_STEPS + 1)} gain_S = {t: [] for t in range(1, NUM_STEPS + 1)} gain_C = {t: [] for t in range(1, NUM_STEPS + 1)} cost_C = {t: [] for t in range(1, NUM_STEPS + 1)} # capital cost (collects) cost_S = {t: [] for t in range(1, NUM_STEPS + 1)} # steel cost (upgrades b,d) # ---- per-city logic ----------------------------------------------- for i in range(N): a_step, a_type = cities[i] # initial type at arrival step init = {"H": isH, "F": isF, "M": isM, "N": isMon} m.Add(init[a_type][i, a_step] == 1) # no upgrades at arrival m.Add(hasA[i, a_step] == 0) m.Add(hasB[i, a_step] == 0) m.Add(hasD[i, a_step] == 0) # vats at arrival m.Add(vE[i, a_step] == 1) m.Add(vB[i, a_step] == 1) m.Add(vS[i, a_step] == 1) # before arrival: everything zero for t in range(1, a_step): for v in (isH, isF, isM, isMon, hasA, hasB, hasD, vE, vB, vS): m.Add(v[i, t] == 0) if t <= NUM_STEPS: for v in ( col, ua, ub, ud, ow, rH, rF, rM, cvE, cvB, cvS, owcvE, owcvB, owcvS, mE, mB, mS, mC, owmE, owmB, owmS, owmC, ): m.Add(v[i, t] == 0) # action + transition logic for active steps for t in range(a_step, NUM_STEPS + 1): P = present[i, t] ren = m.NewBoolVar("") # any renovation this step m.Add(ren == rH[i, t] + rF[i, t] + rM[i, t]) no_ren = m.NewBoolVar("") m.Add(no_ren == 1 - ren) # exactly one action while present m.Add( col[i, t] + ua[i, t] + ub[i, t] + ud[i, t] + ow[i, t] + rH[i, t] + rF[i, t] + rM[i, t] == P ) # renovation must change type (renovate-to-X requires not-X now) m.Add(isH[i, t] == 0).OnlyEnforceIf(rH[i, t]) m.Add(isF[i, t] == 0).OnlyEnforceIf(rF[i, t]) m.Add(isM[i, t] == 0).OnlyEnforceIf(rM[i, t]) # upgrade legality m.Add(hasA[i, t] == 0).OnlyEnforceIf(ua[i, t]) # can't re-acquire m.Add(hasB[i, t] == 0).OnlyEnforceIf(ub[i, t]) m.Add(hasD[i, t] == 0).OnlyEnforceIf(ud[i, t]) m.Add(isF[i, t] == 1).OnlyEnforceIf(ud[i, t]) # d is foundry-only # LLM allowed renovation into metropolis, need to prevent that now m.Add(rM[i, t] == 0) # Monument constraints due to maximizing resources and not wanting enemy to get free upgrades m.Add(col[i, t] == 0).OnlyEnforceIf(isMon[i, t]) m.Add(ua[i, t] == 0).OnlyEnforceIf(isMon[i, t]) m.Add(ub[i, t] == 0).OnlyEnforceIf(isMon[i, t]) m.Add(ud[i, t] == 0).OnlyEnforceIf(isMon[i, t]) m.Add(ow[i, t] == 0).OnlyEnforceIf(isMon[i, t]) # overwork cooldown: city that overworked last step can't collect or overwork if t > a_step: m.Add(col[i, t] == 0).OnlyEnforceIf(ow[i, t - 1]) m.Add(ow[i, t] == 0).OnlyEnforceIf(ow[i, t - 1]) # ---- type transition t -> t+1 ---- m.Add(isH[i, t + 1] == isH[i, t]).OnlyEnforceIf(no_ren) m.Add(isF[i, t + 1] == isF[i, t]).OnlyEnforceIf(no_ren) m.Add(isM[i, t + 1] == isM[i, t]).OnlyEnforceIf(no_ren) for r, target in ((rH[i, t], isH), (rF[i, t], isF), (rM[i, t], isM)): m.Add(target[i, t + 1] == 1).OnlyEnforceIf(r) m.Add(isF[i, t + 1] == 0).OnlyEnforceIf(rH[i, t]) m.Add(isM[i, t + 1] == 0).OnlyEnforceIf(rH[i, t]) m.Add(isH[i, t + 1] == 0).OnlyEnforceIf(rF[i, t]) m.Add(isM[i, t + 1] == 0).OnlyEnforceIf(rF[i, t]) m.Add(isH[i, t + 1] == 0).OnlyEnforceIf(rM[i, t]) m.Add(isF[i, t + 1] == 0).OnlyEnforceIf(rM[i, t]) # ---- upgrade transition t -> t+1 ---- # a, b survive renovation and are monotone m.AddMaxEquality(hasA[i, t + 1], [hasA[i, t], ua[i, t]]) m.AddMaxEquality(hasB[i, t + 1], [hasB[i, t], ub[i, t]]) # d is stripped on renovation, else monotone m.Add(hasD[i, t + 1] == 0).OnlyEnforceIf(ren) d_keep = m.NewBoolVar("") m.AddMaxEquality(d_keep, [hasD[i, t], ud[i, t]]) m.Add(hasD[i, t + 1] == d_keep).OnlyEnforceIf(no_ren) # ---- collect sub-choices ---- fcol = AND(isF[i, t], col[i, t]) mcol = AND(isM[i, t], col[i, t]) hcol = AND(isH[i, t], col[i, t]) # foundry: exactly one vat chosen iff foundry collects m.Add(cvE[i, t] + cvB[i, t] + cvS[i, t] == fcol) # metropolis: pick (2 + hasB) resources (+1 each); else nothing npick = m.NewIntVar(0, 3, "") m.Add(npick == 2 + hasB[i, t]).OnlyEnforceIf(mcol) m.Add(npick == 0).OnlyEnforceIf(mcol.Not()) m.Add(mE[i, t] + mB[i, t] + mS[i, t] + mC[i, t] == npick) # ================= COSTS ================= # foundry & metropolis collect each cost 1 Capital cost_C[t].append(fcol) cost_C[t].append(mcol) # upgrades b,d cost 2 Steel, reduced by 1 if cost-reduction (a) already held for upg in (ub[i, t], ud[i, t]): c = m.NewIntVar(0, 2, "") both = AND(upg, hasA[i, t]) m.Add(c == 2).OnlyEnforceIf(upg, hasA[i, t].Not()) m.Add(c == 1).OnlyEnforceIf(both) m.Add(c == 0).OnlyEnforceIf(upg.Not()) cost_S[t].append(c) # ================= GAINS ================= # Hub collect: +2 Capital (+1 more if collect-bonus b) hub_gain = m.NewIntVar(0, 3, "") m.Add(hub_gain == 2 + hasB[i, t]).OnlyEnforceIf(hcol) m.Add(hub_gain == 0).OnlyEnforceIf(hcol.Not()) gain_C[t].append(hub_gain) # Metropolis collect: +1 per pick gain_E[t].append(mE[i, t]) gain_B[t].append(mB[i, t]) gain_S[t].append(mS[i, t]) gain_C[t].append(mC[i, t]) # Foundry collect: gain chosen vat's value as that resource gEf = m.NewIntVar(0, max_vat, "") gBf = m.NewIntVar(0, max_vat, "") gSf = m.NewIntVar(0, max_vat, "") m.AddMultiplicationEquality(gEf, [cvE[i, t], vE[i, t]]) m.AddMultiplicationEquality(gBf, [cvB[i, t], vB[i, t]]) m.AddMultiplicationEquality(gSf, [cvS[i, t], vS[i, t]]) gain_E[t].append(gEf) gain_B[t].append(gBf) gain_S[t].append(gSf) # Collect Bonus (b): adds +1 to the amount a Collect gives. For a # foundry that means +1 of the collected resource (vat value + 1), # the same uniform "+1 to what Collect gives" rule as Hub/Metro. gain_E[t].append(AND(cvE[i, t], hasB[i, t])) gain_B[t].append(AND(cvB[i, t], hasB[i, t])) gain_S[t].append(AND(cvS[i, t], hasB[i, t])) # ---- overwork sub-choices (double-collect; no Capital cost) ---- fow = AND(isF[i, t], ow[i, t]) mow = AND(isM[i, t], ow[i, t]) how = AND(isH[i, t], ow[i, t]) # foundry overwork: exactly one vat chosen iff foundry overworks m.Add(owcvE[i, t] + owcvB[i, t] + owcvS[i, t] == fow) # metropolis overwork: pick 2*(2 + hasB) resources; else nothing ow_npick = m.NewIntVar(0, 6, "") m.Add(ow_npick == 2 * (2 + hasB[i, t])).OnlyEnforceIf(mow) m.Add(ow_npick == 0).OnlyEnforceIf(mow.Not()) m.Add(owmE[i, t] + owmB[i, t] + owmS[i, t] + owmC[i, t] == ow_npick) # ================= OVERWORK GAINS (2x normal collect) ================= # Hub overwork: 2*(2 + hasB) Capital hub_ow_gain = m.NewIntVar(0, 6, "") m.Add(hub_ow_gain == 2 * (2 + hasB[i, t])).OnlyEnforceIf(how) m.Add(hub_ow_gain == 0).OnlyEnforceIf(how.Not()) gain_C[t].append(hub_ow_gain) # Metropolis overwork: +1 per pick (picks are already doubled via ow_npick) gain_E[t].append(owmE[i, t]) gain_B[t].append(owmB[i, t]) gain_S[t].append(owmS[i, t]) gain_C[t].append(owmC[i, t]) # Foundry overwork: 2 * chosen vat's value owgEf = m.NewIntVar(0, 2 * max_vat, "") owgBf = m.NewIntVar(0, 2 * max_vat, "") owgSf = m.NewIntVar(0, 2 * max_vat, "") _owE = m.NewIntVar(0, max_vat, "") _owB = m.NewIntVar(0, max_vat, "") _owS = m.NewIntVar(0, max_vat, "") m.AddMultiplicationEquality(_owE, [owcvE[i, t], vE[i, t]]) m.AddMultiplicationEquality(_owB, [owcvB[i, t], vB[i, t]]) m.AddMultiplicationEquality(_owS, [owcvS[i, t], vS[i, t]]) m.Add(owgEf == 2 * _owE) m.Add(owgBf == 2 * _owB) m.Add(owgSf == 2 * _owS) gain_E[t].append(owgEf) gain_B[t].append(owgBf) gain_S[t].append(owgSf) # Collect Bonus (b) for foundry overwork: 2 * (+1) = +2 of that resource ow_be = AND(owcvE[i, t], hasB[i, t]) ow_bb = AND(owcvB[i, t], hasB[i, t]) ow_bs = AND(owcvS[i, t], hasB[i, t]) gain_E[t].append(ow_be) gain_E[t].append(ow_be) # added twice == *2 gain_B[t].append(ow_bb) gain_B[t].append(ow_bb) gain_S[t].append(ow_bs) gain_S[t].append(ow_bs) # ---- vat update producing vat[i, t+1] ---- # increment added to the two non-collected vats (1, or 2 with upgrade d) inc = m.NewIntVar(1, 2, "") m.Add(inc == 1 + hasD[i, t]) # vat_next = result of this step's action (only meaningful if foundry) vEn = m.NewIntVar(0, max_vat, "") vBn = m.NewIntVar(0, max_vat, "") vSn = m.NewIntVar(0, max_vat, "") # collect E: E->0, B,S += inc m.Add(vEn == 0).OnlyEnforceIf(cvE[i, t]) m.Add(vBn == vB[i, t] + inc).OnlyEnforceIf(cvE[i, t]) m.Add(vSn == vS[i, t] + inc).OnlyEnforceIf(cvE[i, t]) # collect B m.Add(vBn == 0).OnlyEnforceIf(cvB[i, t]) m.Add(vEn == vE[i, t] + inc).OnlyEnforceIf(cvB[i, t]) m.Add(vSn == vS[i, t] + inc).OnlyEnforceIf(cvB[i, t]) # collect S m.Add(vSn == 0).OnlyEnforceIf(cvS[i, t]) m.Add(vEn == vE[i, t] + inc).OnlyEnforceIf(cvS[i, t]) m.Add(vBn == vB[i, t] + inc).OnlyEnforceIf(cvS[i, t]) # overwork vat transitions: same reset/increment as normal collect m.Add(vEn == 0).OnlyEnforceIf(owcvE[i, t]) m.Add(vBn == vB[i, t] + inc).OnlyEnforceIf(owcvE[i, t]) m.Add(vSn == vS[i, t] + inc).OnlyEnforceIf(owcvE[i, t]) m.Add(vBn == 0).OnlyEnforceIf(owcvB[i, t]) m.Add(vEn == vE[i, t] + inc).OnlyEnforceIf(owcvB[i, t]) m.Add(vSn == vS[i, t] + inc).OnlyEnforceIf(owcvB[i, t]) m.Add(vSn == 0).OnlyEnforceIf(owcvS[i, t]) m.Add(vEn == vE[i, t] + inc).OnlyEnforceIf(owcvS[i, t]) m.Add(vBn == vB[i, t] + inc).OnlyEnforceIf(owcvS[i, t]) # foundry but neither collecting nor overworking: vats unchanged # f_noncollect_noow = isF AND NOT fcol AND NOT fow = isF - fcol - fow f_noncollect_noow = m.NewBoolVar("") m.Add(f_noncollect_noow == isF[i, t] - fcol - fow) m.Add(vEn == vE[i, t]).OnlyEnforceIf(f_noncollect_noow) m.Add(vBn == vB[i, t]).OnlyEnforceIf(f_noncollect_noow) m.Add(vSn == vS[i, t]).OnlyEnforceIf(f_noncollect_noow) # assign vat[i, t+1]: # renovate-to-foundry -> reset to 1 # continuing foundry -> vat_next # otherwise (not foundry next) -> 0 cont_F = AND(isF[i, t], no_ren) # stays a foundry next step for vnext, vn in ( (vE[i, t + 1], vEn), (vB[i, t + 1], vBn), (vS[i, t + 1], vSn), ): m.Add(vnext == 1).OnlyEnforceIf(rF[i, t]) m.Add(vnext == vn).OnlyEnforceIf(cont_F) m.Add(isF[i, t + 1] == 0).OnlyEnforceIf( rF[i, t].Not(), cont_F.Not() ) # tautology guard not_F_next = isF[i, t + 1].Not() m.Add(vE[i, t + 1] == 0).OnlyEnforceIf(not_F_next) m.Add(vB[i, t + 1] == 0).OnlyEnforceIf(not_F_next) m.Add(vS[i, t + 1] == 0).OnlyEnforceIf(not_F_next) # ---- global overwork limit: at most 1 city overworks per step ---- for t in range(1, NUM_STEPS + 1): m.Add(sum(ow[i, t] for i in range(N)) <= 1) # ---- resource pool recursion -------------------------------------- E = {1: m.NewIntVar(initial[0], initial[0], "E1")} B = {1: m.NewIntVar(initial[1], initial[1], "B1")} S = {1: m.NewIntVar(initial[2], initial[2], "S1")} C = {1: m.NewIntVar(initial[3], initial[3], "C1")} for t in range(1, NUM_STEPS + 1): E[t + 1] = m.NewIntVar(0, max_res, f"E_{t + 1}") B[t + 1] = m.NewIntVar(0, max_res, f"B_{t + 1}") S[t + 1] = m.NewIntVar(0, max_res, f"S_{t + 1}") C[t + 1] = m.NewIntVar(0, max_res, f"C_{t + 1}") # costs are paid from the START-of-step pool (must work out before gains) m.Add(C[t] - sum(cost_C[t]) >= 0) m.Add(S[t] - sum(cost_S[t]) >= 0) # next pool = start - costs + gains m.Add(E[t + 1] == E[t] + sum(gain_E[t])) m.Add(B[t + 1] == B[t] + sum(gain_B[t])) m.Add(S[t + 1] == S[t] - sum(cost_S[t]) + sum(gain_S[t])) m.Add(C[t + 1] == C[t] - sum(cost_C[t]) + sum(gain_C[t])) finalE, finalB, finalS = E[NUM_STEPS + 1], B[NUM_STEPS + 1], S[NUM_STEPS + 1] # ---- Phase 1: real per-resource ceilings (fast LINEAR maximisations) ---- # The triple product E*B*S has a very loose relaxation, so proving optimality # directly is slow. We first find the true max each resource can reach on its # own (a linear objective, solved to optimality quickly), then clamp the final # pools to those ceilings. This tightens the product's propagation enough to # prove optimality, and is valid because no feasible solution can exceed an # individual resource's standalone maximum. def _ceiling(var): s = cp_model.CpSolver() s.parameters.max_time_in_seconds = 20.0 s.parameters.num_search_workers = num_workers m.Maximize(var) st = s.Solve(m) return ( int(s.ObjectiveValue()) if st in (cp_model.OPTIMAL, cp_model.FEASIBLE) else max_res ) capE = _ceiling(finalE) capB = _ceiling(finalB) capS = _ceiling(finalS) m.Add(finalE <= capE) m.Add(finalB <= capB) m.Add(finalS <= capS) # ====================================================================== # OBJECTIVE IS SET HERE -- maximise Electrum * Brass * Steel (post step 5) # To change the objective, edit the three "final" pools and/or the product # below. (finalE/finalB/finalS are the pools after step 5's gains.) # ====================================================================== ## NOTE: product can be changed here # prodEB = m.NewIntVar(0, capE * capB, "prodEB") # m.AddMultiplicationEquality(prodEB, [finalE, finalB]) # obj = m.NewIntVar(0, capE * capB * capS, "obj") # m.AddMultiplicationEquality(obj, [prodEB, finalS]) # m.Maximize(obj) # Linear objective instead (it sucks) # m.Maximize(finalE + finalB + finalS) # New Product def Eprod(v): return v * v def Bprod(v): return v def Sprod(v): return v * v prodEE = m.NewIntVar(0, Eprod(capE), "prodEE") m.AddMultiplicationEquality(prodEE, [finalE]) prodSS = m.NewIntVar(0, Sprod(capS), "prodSS") m.AddMultiplicationEquality(prodSS, [finalS]) prodBB = m.NewIntVar(0, Bprod(capB), "prodBB") m.AddMultiplicationEquality(prodBB, [finalB]) prodEB = m.NewIntVar(0, Eprod(capE) * Bprod(capB), "prodEB") m.AddMultiplicationEquality(prodEB, [prodEE, prodBB]) obj = m.NewIntVar(0, Eprod(capE) * Bprod(capB) * Sprod(capS), "obj") m.AddMultiplicationEquality(obj, [prodEB, prodSS]) m.Maximize(obj) # ---- Phase 2: solve the product to optimality ---- solver = cp_model.CpSolver() solver.parameters.max_time_in_seconds = time_limit solver.parameters.num_search_workers = num_workers status = solver.Solve( m, printer.IntermediateSolutionPrinter( {"electrum": finalE, "brass": finalB, "steel": finalS} ), ) if verbose: print(f"(resource ceilings used: E<={capE} B<={capB} S<={capS})") _report( solver, status, cities, N, isH, isF, isM, isMon, col, ua, ub, ud, ow, rH, rF, rM, cvE, cvB, cvS, owcvE, owcvB, owcvS, mE, mB, mS, mC, owmE, owmB, owmS, owmC, hasA, hasB, hasD, vE, vB, vS, E, B, S, C, finalE, finalB, finalS, ) return solver, status def _report( solver, status, cities, N, isH, isF, isM, isMon, col, ua, ub, ud, ow, rH, rF, rM, cvE, cvB, cvS, owcvE, owcvB, owcvS, mE, mB, mS, mC, owmE, owmB, owmS, owmC, hasA, hasB, hasD, vE, vB, vS, E, B, S, C, finalE, finalB, finalS, ): print("status:", solver.StatusName(status)) if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE): return typ_name = {} for i in range(N): for t in range(1, NUM_STEPS + 1): if solver.Value(isH[i, t]): typ_name[i, t] = "Hub" elif solver.Value(isF[i, t]): typ_name[i, t] = "Foundry" elif solver.Value(isM[i, t]): typ_name[i, t] = "Metro" elif solver.Value(isMon[i, t]): typ_name[i, t] = "Monument" else: typ_name[i, t] = "-" def action_str(i, t): if solver.Value(col[i, t]): if typ_name[i, t] == "Foundry": v = ( "E" if solver.Value(cvE[i, t]) else "B" if solver.Value(cvB[i, t]) else "S" ) return f"Collect vat {v}" if typ_name[i, t] == "Metro": picks = [] for nm, var in (("E", mE), ("B", mB), ("S", mS), ("C", mC)): picks += [nm] * solver.Value(var[i, t]) return "Collect {" + ",".join(picks) + "}" return "Collect (+Capital)" if solver.Value(ow[i, t]): if typ_name[i, t] == "Foundry": v = ( "E" if solver.Value(owcvE[i, t]) else "B" if solver.Value(owcvB[i, t]) else "S" ) return f"Overwork vat {v} (2x)" if typ_name[i, t] == "Metro": picks = [] for nm, var in (("E", owmE), ("B", owmB), ("S", owmS), ("C", owmC)): picks += [nm] * solver.Value(var[i, t]) return "Overwork {" + ",".join(picks) + "} (2x)" return "Overwork (+Capital 2x)" if solver.Value(ua[i, t]): return "Upgrade a (cost-reduction)" if solver.Value(ub[i, t]): return "Upgrade b (collect-bonus)" if solver.Value(ud[i, t]): return "Upgrade d (vat-increment)" if solver.Value(rH[i, t]): return "Renovate -> Hub" if solver.Value(rF[i, t]): return "Renovate -> Foundry" if solver.Value(rM[i, t]): return "Renovate -> Metro" return "(inactive)" print("\nPer-step pools (start of step):") print(" step: E B S C") for t in range(1, NUM_STEPS + 2): label = f"after5" if t == NUM_STEPS + 1 else f"start {t}" print( f" {label:>8}: {solver.Value(E[t]):3} {solver.Value(B[t]):3} " f"{solver.Value(S[t]):3} {solver.Value(C[t]):3}" ) print("\nActions:") for i in range(N): a_step, a_type = cities[i] print(f" City {i} (arrives step {a_step} as {a_type}):") for t in range(a_step, NUM_STEPS + 1): ups = "".join( n for n, h in (("a", hasA), ("b", hasB), ("d", hasD)) if solver.Value(h[i, t]) ) extra = ( f" vats(E{solver.Value(vE[i, t])},B{solver.Value(vB[i, t])},S{solver.Value(vS[i, t])})" if typ_name[i, t] == "Foundry" else "" ) print( f" step {t}: [{typ_name[i, t]:7}] {action_str(i, t):28}" f" upg[{ups}]{extra}" ) fe, fb, fs = solver.Value(finalE), solver.Value(finalB), solver.Value(finalS) print( f"\nFINAL E={fe} B={fb} S={fs} product = {fe * fb * fs} sum = {fe + fb + fs}" ) if __name__ == "__main__": solve()