From [1]:
The hour, minute and second hands of this clock are all the same length and move smoothly in a circle. The dial contains hour and minute markers, but the numbers are missing. Therefore, it’s impossible to tell which one of the 12 hour markers belongs to the 12. The two hands on the left are positioned exactly on hour markers, and the hand on the right is positioned between a minute and an hour marker. What time does the clock show?
It is possible to solve this without really any math, but, of course, here I try to model this as a mathematical programming model.
Mathematical Model
| hand | relative position |
|---|---|
| \(h_1\) | \(0\) |
| \(h_2\) | \(2\) |
| \(h_3\) | somewhere inside \((6.8,7)\) |
Notes:
- The symbol \(\mathit{\color{darkred}{relpos}}\) is really a parameter (data), but as I am not exactly sure about \(h_3\), I introduced a feasible interval for that one. The consequence is that this becomes a variable.
- A clock gives rise to modular arithmetic. We have 10:00 + 3 hours = 1:00.
- We can't have a mod operator in a standard mathematical programming model. So the first equation is written as: \[\begin{align}& \mathit{\color{darkred}{hpos}}_{h} = \mathit{\color{darkred}{relpos}}_{h} + \mathit{\color{darkred}{offset}} - 12\color{darkred}\delta_h \\ & \color{darkred}\delta_h \in \{0,1\}\end{align}\]
- We can't really do \(\mathit{\color{darkred}{relpos}}_{h3} \in (6.8,7)\) so we settle for \(\mathit{\color{darkred}{relpos}}_{h3} \in [6.8,7]\).
- The domain of \(\mathit{\color{darkred}{hpos}}_{h}\) is \([0, 11.999]\) to prevent ambiguity at 0 vs 12. There is no hard rule for choosing the tolerance. Too small, and feasibility tolerances in the solver render it useless. Too large: a big hole may make the model infeasible. My choice seems reasonable.
The next task is to model the assignment of hands to their role: indicating hours, minutes and seconds. To form a unique assignment of hands to roles, we can use a fragment of an assignment model: \[\begin{align} &\sum_h \mathit{\color{darkred}{assign}}_{h,r} = 1 && \forall r\\ & \sum_r \mathit{\color{darkred}{assign}}_{h,r} = 1 && \forall h \\ & \mathit{\color{darkred}{assign}}_{h,r} \in \{0,1\}\end{align}\] This gives us a mapping between hands and roles. Feasible levels for this can look like:
---- 86 VARIABLE assign.L assign unique role to hand hour minute second h1 1.000 h2 1.000 h3 1.000
Now we can move from \(\mathit{\color{darkred}{hpos}}_{h}\) to something like \(\mathit{\color{darkred}{rpos}}_{r}\) as follows: \[\mathit{\color{darkred}{rpos}}_{r} = \sum_h \mathit{\color{darkred}{assign}}_{h,r} \cdot \mathit{\color{darkred}{hpos}}_{h}\] Applying this we can see something like:
---- 87 VARIABLE hpos.L final position of hands (relpos + offset mod 12, in decimal hours) h1 10.000, h3 4.833 ---- 87 VARIABLE assign.L assign unique role to hand hour minute second h1 1.000 h2 1.000 h3 1.000 ---- 87 VARIABLE rpos.L position of hr/min/sec hand ∈[0,11.99] hour 4.833, minute 10.000
Note that hand \(h_2\) corresponding to seconds is at position zero (zeros are not printed here).
There is a slight issue here. As we multiply two variables, this formulation would make the model a nonconvex MIQCP (Mixed-Integer Quadratically Constrained). However, we can try a linearization. The idea is that we actually can linearize the product between a binary variable and a (bounded) continuous variable. Consider the case \(\color{darkred}z = \color{darkred}\delta\cdot \color{darkred}x\) with \(\color{darkred}\delta \in \{0,1\}\) (binary) and \(\color{darkred}x \in [0,\color{darkblue}U]\) (continuous), then we can form the set of inequalities: \[\boxed{\begin{aligned} &\color{darkred}z = \color{darkred}\delta\cdot \color{darkred}x\\ & \color{darkred}\delta \in \{0,1\}\\& \color{darkred}x \in [0,\color{darkblue}U] \end{aligned}} \iff \boxed{\begin{aligned} & \color{darkred}z \le \color{darkblue}U \cdot \color{darkred}\delta \\& \color{darkred}z \le \color{darkred}x \\ & \color{darkred}z \ge \color{darkred}x - \color{darkblue}U\cdot (1-\color{darkred}\delta)\\ & \color{darkred}z \in [0,\color{darkblue}U]\\ & \color{darkred}\delta \in \{0,1\}\\& \color{darkred}x \in [0,\color{darkblue}U] \end{aligned}}\] The formulation on the right is completely linear. You can easily verify the correctness of this by substituting \(\color{darkred}\delta=0\) and \(\color{darkred}\delta=1\). It is noted that in our problem we have good, tight bounds on the continuous variables that appear in the multiplication (they are \(\mathit{\color{darkred}{hpos}}_{h} \in [0,11.999]\)). So this linearization technique is perfectly fine to use. Without good bounds, this method can be dangerous. And just to confirm: this is an exact linearization. It is not an approximation. Here is a good example where MINLP formulations can lead to more compact models.
To create proper values for hours, minutes and seconds we can add: \[\begin{align}&\color{darkred}v_{\mathit{second}} = \frac{60}{12}\mathit{\color{darkred}{rpos}}_{\mathit{second}} \\ & \color{darkred}v_{\mathit{minute}}=\frac{60}{12}\mathit{\color{darkred}{rpos}}_{\mathit{minute}}- \frac{1}{60}\color{darkred}v_{\mathit{second}} \\ & \color{darkred}v_{\mathit{hour}} = \mathit{\color{darkred}{rpos}}_{\mathit{hour}} - \frac{1}{60}\color{darkred}v_{\mathit{minute}} - \frac{1}{60\cdot 60}\color{darkred}v_{\mathit{second}} \\ & \color{darkred}v_{\mathit{hour}} \in \{0,1,\dots,11\} \\ & \color{darkred}v_{\mathit{minute}},\color{darkred}v_{\mathit{second}} \in \{0,1,\dots,59\}\end{align}\] This also establishes relationships between hours, minutes and seconds, i.e. we get a well-defined representation of time.
We don't have an objective: this is a feasibility problem.
With this, we have solved the problem:
---- 95 VARIABLE v.L (integer) value in hours,minutes and seconds hour 4.000, minute 50.000
so our time would be 4:50.
We started this discussion using imprecise data for the third hand. That may make you wonder: does this allow for alternative solutions? A possible way to assess this, is to form two optimization problems: one minimizing the "time" (\(\mathit{\color{darkred}{rpos}}_{\mathit{hour}}\)) and the other one maximizing this. This is essentially a form of sensitivity analysis. This exercise did not show any difference.
Although this is actually more work than directly solving the problem by hand, it is still useful. We have discussed several interesting modeling issues (modular arithmetic, linearization of a product). And also, developing a model forces you to have a look at the problem a bit differently, and a bit more precisely. Sometimes, a model is not that important, but the modeling process is still a good way to gain a better understanding of the problem.
GAMS Models (nonlinear and linear)
$onText
Clock problem
see:
https://www.scientificamerican.com/game/math-puzzle-find-the-time/
$offText
Set
h 'hand' /h1,h2,h3/
r 'role' /hour,minute,second/
;
variable
* hands h1,h2,h3
relpos(h) 'relative position of hands (define h1=0)'
offset 'from relpos'
hpos(h) 'final position of hands (relpos + offset mod 12, in decimal hours)'
modulo(h) 'binary variable to implement mod 12'
* role hour,minute,second
assign(h,r) 'assign unique role to hand'
rpos(r) 'position of hr/min/sec hand ∈[0,11.99]'
* proper clock values
v(r) '(integer) value in hours,minutes and seconds'
* objective
z 'dummy objective variable'
;
integer variable offset,v;
binary variable modulo, assign;
offset.up = 11;
hpos.lo(h) = 0;
hpos.up(h) = 11.999;
v.lo(r) = 0;
v.up(r) = 59;
v.up('hour') = 11;
* read relative positions from the clock:
* h1 = 0, h2 = 2, h3 ∈ [6.8,7]
table relposdata(h,*) 'relative position (interval) in units of hours (define h1=0)'
lo up
h1 0 0
h2 2 2
h3 6.8 7
;
relpos.lo(h) = relposdata(h,'lo');
relpos.up(h) = relposdata(h,'up');
equations
obj 'dummy objective'
e_hpos(h) 'position of hands (in decimal hours)'
e_assign1 'unique assignment hand <=> role'
e_assign2 'unique assignment hand <=> role'
assignpos 'calc rpos from hpos (nonlinear)'
* unit conversions
seconds 'calc seconds from rpos'
minutes 'calc minutes from rpos'
hours 'calc hours from rpos'
;
* label(h) = (relpos(h) + offset) mod 12
e_hpos(h).. hpos(h) =e= relpos(h) + offset - modulo(h)*12;
* unique assignment hand <=> role
e_assign1(r).. sum(h,assign(h,r)) =e= 1;
e_assign2(h).. sum(r,assign(h,r)) =e= 1;
assignpos(r).. rpos(r) =e= sum(h, assign(h,r)*hpos(h));
* proper clock values
seconds.. v('second') =e= rpos('second')*60/12;
minutes.. v('minute') =e= rpos('minute')*60/12 - v('second')/60;
hours.. v('hour') =e= rpos('hour') - v('minute')/60 - v('second')/60/60;
obj.. z =e= 0;
*------------------------------------------------
* nonconvex quadratic model
*------------------------------------------------
model m1 /all/;
option miqcp=baron;
solve m1 minimizing z using miqcp;
display relpos.l,offset.l,hpos.l,assign.l,rpos.l,v.l;
*------------------------------------------------
* linearized model
*------------------------------------------------
positive variables prd(h,r) 'product prd(h,r)=assign(h,r)*hpos(h)';
prd.up(h,r) = hpos.up(h);
equations
prod1(h,r) 'linearization of product'
prod2(h,r) 'linearization of product'
prod3(h,r) 'linearization of product'
assignpos2(r) 'linearization of product'
;
prod1(h,r).. prd(h,r) =l= hpos.up(h)*assign(h,r);
prod2(h,r).. prd(h,r) =l= hpos(h);
prod3(h,r).. prd(h,r) =g= hpos(h)-hpos.up(h)*(1-assign(h,r));
assignpos2(r).. rpos(r) =e= sum(h, prd(h,r));
model m2 /m1-assignpos,prod1,prod2,prod3,assignpos2/;
solve m2 minimizing z using mip;
display relpos.l,offset.l,hpos.l,assign.l,rpos.l,v.l;
ChatGPT
It is interesting to see what the popular LLMs produce when asked to formulate a Mixed-Integer Programming problem for this problem. I get all very impressive-looking models, but most of them were incorrect. The best was:
Pyomo Model generated by ChatGPT
import pyomo.environ as pyo # ----------------------------- # Approximated observed angles # ----------------------------- a = {1: 300.0, # minute hand 2: 0.0, # second hand 3: 145.0} # hour hand # Hands known to be on hour markers (from picture) L = [1, 2] R = [3] # Big-M constant Mbig = 1000 model = pyo.ConcreteModel() # ----------------------------------------- # Sets # ----------------------------------------- model.I = pyo.Set(initialize=[1,2,3]) # observed hands model.K = pyo.Set(initialize=['h','m','s']) # roles model.L = pyo.Set(initialize=L) model.R = pyo.Set(initialize=R) # ----------------------------------------- # Decision variables # ----------------------------------------- # Time variables model.H = pyo.Var(domain=pyo.Integers, bounds=(0,11)) model.M = pyo.Var(domain=pyo.Integers, bounds=(0,59)) model.S = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0,60)) # Dial rotation (which hour is "12") model.r = pyo.Var(domain=pyo.Integers, bounds=(0,11)) # Assignment binaries model.y = pyo.Var(model.I, model.K, domain=pyo.Binary) # Angle variables model.theta = pyo.Var(model.K, bounds=(0,360)) # Modulo-integer variables (at most ±360) model.z = pyo.Var(model.I, model.K, domain=pyo.Integers, bounds=(-1,1)) # Hour-marker position for left-side hands model.hpos = pyo.Var(model.L, domain=pyo.Integers, bounds=(0,11)) model.q = pyo.Var(model.L, domain=pyo.Integers, bounds=(-1,1)) # ----------------------------------------- # Clock kinematic constraints # ----------------------------------------- def hour_angle(m): return m.theta['h'] == 30*m.H + 0.5*m.M + (1/120)*m.S model.c_hour = pyo.Constraint(rule=hour_angle) def minute_angle(m): return m.theta['m'] == 6*m.M + 0.1*m.S model.c_minute = pyo.Constraint(rule=minute_angle) def second_angle(m): return m.theta['s'] == 6*m.S model.c_second = pyo.Constraint(rule=second_angle) # ----------------------------------------- # Assignment (permutation) constraints # ----------------------------------------- def one_role_each_hand(m,i): return sum(m.y[i,k] for k in m.K) == 1 model.c_assign_hand = pyo.Constraint(model.I, rule=one_role_each_hand) def one_hand_each_role(m,k): return sum(m.y[i,k] for i in m.I) == 1 model.c_assign_role = pyo.Constraint(model.K, rule=one_hand_each_role) # ----------------------------------------- # Angle matching constraints with modulo # ----------------------------------------- def angle_match(m,i,k): # (a_i - 30*r) - theta_k + 360*z = 0 if assigned return (a[i] - 30*m.r) - m.theta[k] + 360*m.z[i,k] <= Mbig*(1 - m.y[i,k]) model.c_angle_upper = pyo.Constraint(model.I, model.K, rule=angle_match) def angle_match_lower(m,i,k): return (a[i] - 30*m.r) - m.theta[k] + 360*m.z[i,k] >= -Mbig*(1 - m.y[i,k]) model.c_angle_lower = pyo.Constraint(model.I, model.K, rule=angle_match_lower) # ----------------------------------------- # Hands exactly on hour markers (left side) # ----------------------------------------- def on_hour_marker(m, i): # (a_i - 30*r) = 30*hpos_i - 360*q_i return (a[i] - 30*m.r) - 30*m.hpos[i] + 360*m.q[i] == 0 model.c_marker = pyo.Constraint(model.L, rule=on_hour_marker) # ----------------------------------------- # Hand not on marker (right side) # (simple version: forbid exact equality) # ----------------------------------------- epsilon = 1e-3 def off_marker(m, i): # |(a_i - 30*r) mod 30| >= epsilon # Use linear disjunction: forbid exact multiples of 30° return abs(a[i] - 30*m.r - 30*round((a[i] - 30*m.r)/30)) >= epsilon # In Pyomo we must skip nonlinear abs — so instead we simply *omit* a marker constraint # for i in R. Solver will place this hand at non-marker automatically. # (No explicit constraint needed.) # ----------------------------------------- # Objective: feasibility only # ----------------------------------------- model.obj = pyo.Objective(expr=0) # ----------------------------------------- # Solve # ----------------------------------------- solver = pyo.SolverFactory('cbc') # choose your solver res = solver.solve(model, tee=True) print("Hour:", pyo.value(model.H)) print("Minute:", pyo.value(model.M)) print("Second:", pyo.value(model.S)) print("Rotation r:", pyo.value(model.r)) print("Assignments:", { (i,k): pyo.value(model.y[i,k]) for i in model.I for k in model.K })
ChatGPT says this will produce:
which is close. We can fix this by repairing the range for the variable representing the seconds (from \([0,60]\) to \([0,59.9]\)).
- I work with relative positions (expressed in "hours"). ChatGPT works with angles. This is more of a unit issue than anything else.
- My linearization is more streamlined.
- Fewer constraints.
- ChatGPT uses a big M formulation with \(M=1,000\). I am more careful and use tighter bounds.
- ChatGPT uses integers for hours and minutes, but a continuous variable for seconds. I use integers for all.
- I am using an explicit interval to handle the slightly uncertain value of the third hand. ChatGPT is using an estimate for the angles.
My standard objection against copy-pasting code from LLMs is that actually putzing around with a model (including making mistakes) is very useful to learn about the problem, the model, and its behavior. Missing out on this is really a disadvantage.
References
- Math Puzzle: Find the Time, https://www.scientificamerican.com/game/math-puzzle-find-the-time/
- Original from: Welche Uhrzeit ist gesucht?, https://www.spektrum.de/raetsel/welche-uhrzeit-ist-gesucht/2266965
No comments:
Post a Comment