Here is the problem:
I tried Python's pyomo library but can not find a solution as there is no modulus operation to decide if the number can be divisible by 13.
This can be modeled fairly straight-forwardly using the MiniZinc modelling language in a similar way to the model done by Erwin Kalvelagen in the comments. The main idea is to build up the factors for each number, and to multiply it all. MiniZinc has native support for the mod operator, although not all solvers can handle it.
There are 108 solutions to your instance found using the built-in Gecode 6.1.1 solver in the MiniZinc IDE version 2.4.3. Neither the Chuffed nor the COIN-BC built-in solvers could handle the model.
The main drawback of this way to model it, is that it requires building up the full value of the concatenation. Since solvers typically have a limited range of values that can be represented, this puts an upper limit on the number and sizes of values. For example, Gecode can handle essentially 32 bit integers as values.
The full model for your problem is given here
include "globals.mzn";
%
% Problem data
%
% The set of allowed values
set of int: Values = {1, 10, 20, 40, 100, 200};
% The number of value to concatenate
int: n = 4;
%
% Derived data
%
% The available values as an index set
set of int: Index = 1..card(Values);
% The values as an array
array[Index] of int: values = [v | v in Values];
% The length of values in decimal notation
array[Index] of int: lengths = [ceil(log10(v+1)) | v in Values];
% A mapping table from values to the widths of the numbers
array[Index, 1..2] of int: value_to_widths = array2d(Index, 1..2, [
[values[i], 10^lengths[i]][v_or_w]
| i in Index, v_or_w in 1..2]);
% The numbers as an index set
set of int: N = 1..4;
%
% Variables
%
% The numbers to choose
array[N] of var Values: numbers;
% The widths of the numbers
array[N] of var int: widths;
% The factors to use for concatenation
array[1..n] of var int: factors;
% The concatenation of the numbers
var int: concatenation;
%
% Constraints
%
% For each number, find the widths of it using the mapping table
constraint forall (i in N) (
table([numbers[i], widths[i]], value_to_widths)
);
% Find the factors to multiply numbers with for the concatenation
constraint factors[n] == 1;
constraint forall (i in 1..n-1) (factors[i] == widths[i+1] * factors[i+1]);
% Concatenate the numbers using the factors
constraint concatenation == sum(i in N) (numbers[i] * factors[i]);
% Problem constraint, the concatenation should be evenly divisible by 13
constraint concatenation mod 13 == 0;
%
% Solve and output
%
% Find solutions using standard search
% 636 failures to find all solutions using Gecode 6.1.1
%solve satisfy;
% Find solution by setting the least significant number first
% 121 failures to find all solutions using Gecode 6.1.1
solve :: int_search(reverse(numbers), input_order, indomain_min) satisfy;
% Find the smallest concatenation using the numbers
%solve minimize concatenation;
% Find the largest concatenation using the numbers
%solve maximize concatenation;
% Output both data and variables
output [
"values=", show(values), "\n",
"lengths=", show(lengths), "\n",
"v_to_w={ "] ++ [show(value_to_widths[i, 1]) ++ "->" ++ show(value_to_widths[i, 2]) ++ " "|i in Index] ++ ["}\n",
"numbers=", show(numbers), "\n",
"widths=", show(widths), "\n",
"factors=", show(factors), "\n",
"concatenation=", show(concatenation), "\n",
];
With the following output including the first and last solution when finding all solutions, and printing statistics
values=[1, 10, 20, 40, 100, 200]
lengths=[1, 2, 2, 2, 3, 3]
v_to_w={ 1->10 10->100 20->100 40->100 100->1000 200->1000 }
numbers=[20, 1, 1, 1]
widths=[100, 10, 10, 10]
factors=[1000, 100, 10, 1]
concatenation=20111
----------
[... lots of more solutions printed ...]
----------
values=[1, 10, 20, 40, 100, 200]
lengths=[1, 2, 2, 2, 3, 3]
v_to_w={ 1->10 10->100 20->100 40->100 100->1000 200->1000 }
numbers=[10, 40, 200, 200]
widths=[100, 100, 1000, 1000]
factors=[100000000, 1000000, 1000, 1]
concatenation=1040200200
----------
==========
%%%mzn-stat: initTime=0.000603
%%%mzn-stat: solveTime=0.003856
%%%mzn-stat: solutions=108
%%%mzn-stat: variables=15
%%%mzn-stat: propagators=12
%%%mzn-stat: propagations=8098
%%%mzn-stat: nodes=457
%%%mzn-stat: failures=121
%%%mzn-stat: restarts=0
%%%mzn-stat: peakDepth=7
%%%mzn-stat-end
Finished in 253msec