File: ef230-2023-08/www/examples/sudoku_solver2.m Download
function A=sudoku_solver2(P)
% This is a version of the example provided in the documention
% that shows how to solve a a Sudoku puzzle using MATLAB's new optimization
% command intlinprog introduced in R2014a
% May 8, 2014 W. Schleter EF 230 University of Tennessee
%
% input should be a 9,9 matrix of the puzzle with zeros representing empty boxes
% output will be the solved puzzle
close all;
format compact;

if ~exist('intlinprog')
    error('This example requires release R2014a or later with the optimization toolbox');
end

if nargin==0
    P = getDefaultPuzzle();
end

drawGrid();
drawNums(P,'k');

lb = zeros(9,9,9);
ub = lb+1;
numvars = 9^3; % number of variables;
numcons = 4*(9^2); % number of constraints i.e. a single value in each row, col, dep, and square
Aeq = zeros(numcons,numvars);
beq = ones(numcons,1);

icon = 0;
for i=1:9
    for j=1:9
        T = lb;
        T(:,i,j) = 1;
        icon=icon+1;
        Aeq(icon,:) = T(:);
        
        T = lb;
        T(i,:,j) = 1;
        icon=icon+1;
        Aeq(icon,:) = T(:);

        T = lb;
        T(i,j,:) = 1;
        icon=icon+1;
        Aeq(icon,:) = T(:);        
    end
end

for i=0:3:6
    for j=0:3:6
        for k=1:9
            T = lb;
            T(i+(1:3),j+(1:3),k) = 1;
            icon=icon+1;
            Aeq(icon,:) = T(:);
        end
    end
end

f=1:numvars; % doesn't matter what f is, just needs to be defined
% f=randi([0 1],numvars,1); this also works
% f=zeros(numvars,1); this also works

ic=1:numvars; % indices of the unknowns that should be constrained to integer values

g=find(P>0); % the indices of the given numbers in the puzzle
lb(g+81*(P(g)-1))=1; % set the lower bound to 1 for the given numbers
% the loop below is equivalent to the line above
% for i=1:length(g)
%     i2=g(i);
%     [r,c] =ind2sub([9,9],i2);
%     lb(r,c,P(i2))=1
% end

tic
[AA,fval,exitflag,output]=intlinprog(f,ic,[],[],Aeq,beq,lb,ub);
toc
A=zeros(9,9);
for i=(find(AA>0))'
    [r,c,d] = ind2sub([9,9,9],i);
    A(r,c)=d;
end
drawNums(A,'b');
drawNums(P,'k');
end % sudoku_solver


function drawGrid()
clf
%axis off;
axis([-1 10 -1 10]);
axis equal
set(gca,'ydir','reverse');
hold all
w = [3 1 1 3 1 1 3 1 1 3]; %line widths
for i=0:9
    plot([i i],[0 9],'g-','linewidth',w(i+1)); % vertical
    plot([0 9],[i i],'g-','linewidth',w(i+1)); % horizontal
end
drawnow
end % drawGrid

function drawNums(S,color)
for r=1:9
    for c=1:9
        if S(r,c)>0
            text(c-0.5,r-0.5,num2str(S(r,c),1),'color',color,'horizontal','center','vertical','middle','fontsize',18);
        end
    end
end
drawnow
end % drawNums

function S = getDefaultPuzzle()
S = [ 0 2 0 0 3 0 0 4 0
      6 0 0 0 0 0 0 0 3
      0 0 4 0 0 0 5 0 0
      0 0 0 8 0 6 0 0 0
      8 0 0 0 1 0 0 0 6
      0 0 0 7 0 5 0 0 0
      0 0 7 0 0 0 6 0 0
      4 0 0 0 0 0 0 0 8
      0 3 0 0 4 0 0 2 0 ];
end

function S = getDefaultPuzzle2()
S = [ 0 0 0 6 0 9 0 0 0
      0 3 0 0 4 2 8 0 0
      0 9 7 0 0 0 0 5 0
      0 5 8 2 0 0 0 0 0
      0 0 6 0 0 0 7 0 0
      0 0 0 0 0 8 4 2 0
      0 8 0 0 0 0 3 1 0
      0 0 1 9 3 0 0 4 0
      0 0 0 4 0 1 0 0 0]
end % getDefaultPuzzle